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Abstract 

The purpose of this paper is to give an overview in the realm of numerical computations of 
polydispersed turbulent two-phase flows, using a mean- field/PDF approach. In this approach, 
the numerical solution is obtained by resorting to a hybrid method where the mean fluid 
properties are computed by solving mean-field (RANS) equations with a classical finite volume 
procedure whereas the local instantaneous properties of the particles are determined by solving 
stochastic differential equations (SDEs). The fundamentals of the general formalism are 
recalled and particular attention is focused on a specific theoretical issue: the treatment of 
the multiscale character of the dynamics of the discrete particles, that is the consistency of 
the system of SDEs in asymptotic cases. Then, the main lines of the particle/mesh algorithm 
are given and some specific problems, related to the integration of the SDEs, are discussed, 
for example, issues related to the specificity of the treatment of the averaging and projection 
operators, the time integration of the SDEs (weak numerical schemes consistent with all 
asymptotic cases), and the computation of the source terms. Practical simulations, for three 
different flows, are performed in order to demonstrate the ability of both the models and the 
numerics to cope with the stringent specificities of polydispersed turbulent two-phase flows. 
Key words: dispersed two-phase flows, turbulence, PDF methods, particle/mesh algorithms. 
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1 Introduction 



Two-phase flows are relatively easy to observe: to get a first inkling, one can think of throwing 
small light particles (which then play the role of tracer particles) into a turbulent flow such as a 
rapid river or a plume coming out of a chimney. The small solid particles reveal the intricate and 
complex features of turbulent flows: understanding and modelling these features, i.e. single-phase 
turbulent flow modelling, is the subject of extensive research pQ. If one introduces larger and larger 
particles in the flow, more complex phenomena take place: the behaviour of the particles will reflect 
the interplay between the main physical mechanisms, such as particle inertia and turbulence of the 
carrying flow. Then, eventually, when particles become large-enough, the effect of the fluid may 
become negligible with respect to particle inertia. Thus, turbulent fluid-particle flow modelling 
appears as a link between subjects such as turbulence and granular flows [2]. In general, two-phase 
flows are even more complex since, in the case of air and water for example, different configurations 
of the interface between the two phases may be present. Yet, in the present study, attention will 
be focused on the motion of particles embedded in a turbulent fluid, i.e. polydispersed turbulent 
two-phase flows, where the geometrical configuration does not change. 

Polydispersed turbulent two-phase flows are found in numerous environmental and industrial 
processes, very often in contexts that involve additional issues, for example chemical and com- 
bustion ones. Therefore, modelling these flows raises very difficult theoretical questions and, at 
the same time, one has to provide answers to what we can refer to as engineering concerns. As 
a consequence, a theoretical and numerical model represents an attempt to find a satisfactory 
compromise between these sometimes conflicting expectations. Before trying to outline what is 
meant by satisfactory, let us describe the characteristics of the polydispersed turbulent two-phase 
flows we consider here. 

In the present study, only non-reacting incompressible fluid-particle flows, with no collisions 
between particles, are investigated (particle dispersion and turbulence modulation induced by the 
presence of the particles are the physical mechanisms under consideration). This is not a strict 
limitation of the approach that will be adopted since, as mentioned below, the probability density 
function (PDF) formalism that shall be followed is precisely well-suited for the extension to more 
complex physics, such as combustion. However, for the sake of simplicity, we limit ourselves to 
the core physics embodied by particle dynamics. In addition, only the case of solid heavy particles 
is treated, i.e. the density of the particles is much greater than that of the fluid, p p ^> pf. This 
hypothesis simplifies the equation of motion of the discrete particles in a turbulent flow which, 
retaining only drag and gravity forces, can be written as: 



In these equations, U s (f) = U(t,x p (f)) is the fluid velocity "seen", i.e. the fluid velocity sampled 
along the particle trajectory x p (£), where U(i, x) is the local instantaneous (Eulcrian) fluid velocity 
field. The particle relaxation time, t p , is defined as 



where the local instantaneous relative velocity is U r (t) = U p (i) — U s (i). The drag coefficient, Cd, 
is a non-linear function of the particle-based Reynolds number, Re p = d p \\J r \/vf, which means 
that Cd is a non-linear function of the particle diameter, d p , [3]. This last point represents a 
major theoretical difficulty for a statistical treatment since we do not consider mono-dispersed 
two-phase flows (where d p is constant), but polydispersed two-phase flows where d p covers a range 
of possible values (from very light particles acting as fluid tracers to high-inertia particles in the 
ballistic regime, where the effect of the fluid on the particle dynamics can be neglected). In the 
particle dynamical equations, it is important to note that we are dealing with the instantaneous 
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fluid velocities, U(£,x p (£)). Yet, for high-Reynolds turbulent flows, which are the most common 
ones, such an information is not available due to the very large number of degrees of freedom of the 
turbulent flows [3]. A modelling step is necessary and most models adopt a statistical approach 
where only some limited information is sought for the fluid fields whereas particles are tracked 
individually. In practical models, this information consists in, for the fluid, the first two velocity 
moments, as in Rij — e models [5], or even filtered velocity fields as in LES calculations, e.g. 
In the present work, — e models (RANS equations) will be used in practical computations, but 
the PDF approach (for the particles) to come is fully compatible with other approaches for the 
fluid, for example LES. 

As indicated above, in order to track the particles, a satisfactory model must be built for the 
evaluation of the particle properties, cf. Eqs. ([T}. By satisfactory, it is meant a model which has 
the following properties: 

(i) the model treats the important phenomena, such as convection and the polydispersed nature 
of the particles, without approximation, 

(ii) the approach is naturally set into a general formalism which allows additional variables, for 
more complex physics such as combustion issues, to be directly introduced, 

(iii) the complete theoretical model must be tractable in complex geometries and applicable for 
engineering problems. 

The first two issues have been addressed in a previous review work [8] where a PDF approach 
has been developed. In practice, the PDF approach has the form of a particle stochastic method 
where the velocity of the fluid seen, U s (£), is modelled as a stochastic diffusion process, i.e. the 
dynamics of the particles are calculated from SDEs (Stochastic Differential Equations), the so- 
called Langevin equations [5]. In polydispersed two-phase flows, a particle point of view seems 
rather natural, given the physics considered. Yet, the particles which are to be simulated represent 
samples of the pdf and should not be confused with real particles. Within the PDF formalism, this 
particle point of view is helpful to build the theoretical model and, at the same time, represents 
directly a discrete formulation of the model. However, in order to devise a consistent framework, 
it is important to separate the two steps by formulating the model in continuous time before 
addressing the questions of numerical methods for practical computations. 

The purpose of the present review is to address point (iii) above, and to discuss the general 
numerical methodology used for particle stochastic or PDF models for polydispersed turbulent 
two-phase flows. More specifically, it is aimed at providing answers to several interrogations: 

(a) what do the stochastic particles represent? 

(b) how do we compute the stochastic differential equations? 

(c) what are the various difficulties and sources of numerical errors in the complete numerical 
method? 

More than trying to present definitive answers to the questions of what numerical scheme should 
be used, the objective is to propose a general numerical approach and to show how PDF models, 
Langevin stochastic equations, particle/mesh and dynamical Monte Carlo methods are closely 
connected and actually represent different translations of the same idea. Within that context, a 
major goal is to emphasise that, although numerical schemes are separated from the construction 
of the theoretical model, they cannot be addressed only from a mathematical point of view. Indeed, 
it is important that they reflect the physical properties of the continuous stochastic model, namely 
the multiscale character of the Langevin equations presented in Section [3] 

The paper is organised as follows. The mathematical background on PDF equations and 
stochastic diffusion processes is recalled in Section[2] General and state-of-the-art Langevin models 
for polydispersed turbulent two-phase flow modelling are briefly presented in Section [3] A central 
point is the analysis of the multiscale properties of the Langevin equations, and the expression 
of the various physical limits when characteristic timescales become negligible with respect to the 
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observation timescale, cf. Section @] This analysis serves as a guideline for the development of the 
numerical model in Section[Sl that contains both particle/mesh and time-integration issues. Then, 
a discussion is given on specific issues related to two-way coupling, Section |BJ Several numerical 
applications representative of practical concerns are proposed in Section [7J 



2 General formalism 

The general formalism on which the derivation of the system of equations (RANS equations for 
the fluid and SDEs for the discrete particles) relies, is based on the Lagrangian point of view: the 
system (the fluid-particle mixture) is treated as an ensemble of fluid and discrete particles. The 
discretisation of a continuous medium (the fluid) with particles is not a natural step, but it is a 
practical way, in the frame of the probabilistic formalism briefly outlined here, to treat important 
physical phenomena without approximation 8 . In the present formalism, a fluid particle is an 
independent sample of the flow, with a given pdf. Physically, a fluid particle can be seen as a 
small element of fluid whose characteristic length scale is much larger than the molecular mean 
free path and much smaller than the Kolmogorov length scale. The fluid particles have a mass 
to/, a volume Vt and a velocity that equals the fluid velocity field at the location of the particle, 
U / (t) = U(t,x / (t)). 

2.1 Statistical approach 

Let us consider an ensemble composed by Nf fluid particles and N p discrete particles interacting 
through forces that can be expressed as functions, or functionals, of the variables attached to each 
particle, e.g., I variables for the fluid particles and q variables for the discrete particles (these 
variables can be, for example, position, velocity, . . . ). All available information is then contained 
in the following state vector: 

Z(t) = {Z) A {t), Zl(t) ; . . . ; Z#(t), . . . , Z% (*); Z^(t), ... 

...,Z^ q (t);...;Z^(t),...,Z^(t)}, 

where Zl . (t) represents the variable j attached to the fluid particle labelled i and Z^ j (t) represents 
the variable j attached to the discrete particle labelled i. The dimension of the state vector is then 
INf + qN p . Let us suppose that the dynamical behaviour of the closed system can be described 
in terms of ordinary differential equations (the Navier-Stokes equations, in Lagrangian form, for 
the fluid particles and for the discrete particles, the equation of motion of a single particle in a 
turbulent fluid-particle mixture), i.e. 

^ = A(t,Z W ). (4) 

Here, it is assumed that, in the Navier-Stokes equations, the local instantaneous pressure gradient, 
the viscous forces and the source term (due to the force exerted by the discrete particles on the 
fluid) can be expressed as functionals of the state vector Z(t). In sample space, this system of 
ODEs (Ordinary Differential Equations) corresponds to the Liouville equation [9] 

^M + ^( A(t)Z ) Kt;z )) =0 , (5) 

where p(t;z), the associated pdf, represents the probability to observe at time t the system in 
state z. In the present paper, we distinguish between physical space, Z, and sample space, z. A 
distinction is also made, for the pdf, between parameters and variables by separating them with 
a semi-colon, i.e. (t;z). 

In practice, the number of degrees of freedom of such a system is huge (turbulent flow with 
a large number of particles) and one has to resort to a contracted description in order to come 



5 



up with a model that can be simulated with modern computer technology. For single-phase 
turbulent reactive flows, a one-point pdf, p(t;z^), is often retained [101 HI]- F° r the description 
of the dynamics of the discrete particles in turbulent dispersed two-phase flows, a one-point pdf, 
pit; z J p ), is also encountered [SJIH1H2]- In this work, as we shall see in Section I3~2l both approaches 
are gathered in the form of a two-point pdf, p{t; z l ^,z p ) and the associated reduced state vector 
(henceforth denoted by superscript r) is 

Z r (t) ={Z ftl (t),..., Z ftl (t), Z pA (t),..., Z p . q (t)}. (6) 

The time evolution equations, in physical space, for this sub-system have the form 

dZ r (t) 



dt 



= A(t,Z r (t),Y(t)), (7) 



where there is a dependence on the external variable Y(t) (related to the particles not contained 
in Z r {t) as only pairs of particles, a fluid one and a discrete one, are under consideration). In 
sample space, the marginal pdf p r (t; z r ) verifies 

+ A [(A |z'-)^;z'-)]=0 ! (8) 
where the conditional expectation is given by 

(A|z'">= [ A(t,z r ,y)p(t;y\z r )dy 

Eq. (|8|) is now unclosed, showing that a reduced description of a system implies a loss of informa- 
tion and thus the necessity to introduce a model. 

A practical way to close the system is to resort to stochastic differential equations (SDEs), as it 
shall be briefly explained in Section 13.11 Further detailed explanations for this move can be found 
in Refs. [5] and [13]. The stochastic differential equations treated in this work have the following 
form (Z r (t) is called a diffusion process) 

dZ\{t) = Ai(t, Z r (t)) dt + Bij(t, Z r {t)) dWj(t), (10) 

where W(t) = (W±(t), . . . , Wd(t)) is a set of independent Wiener processes [14] and d = l + q is the 
dimension of the reduced state vector. These equations are often referred to as Langevin equations 
in the physical literature [5]. In Eq. (fTU|). A — (Ai) is called the drift vector and B = (.By) is 
the diffusion matrix. SDEs require a strict mathematical definition of the stochastic integral as 
it shall be explained in Section 12.21 If one adopts the Ito definition of the stochastic integral (see 
Section [2~2l in Eq. (TlOl) . it can be shown, see [9] for example, that the corresponding equation in 
sample space for p r (t; z r ) is the Fokker-Planck equation 

Brf(t-f r \ 8 1 r) 2 

^H = "^ [ ^ ( '' zr)pr(i;zr)] + 2 9<a4 [ ^ (< ' zr)pr( * ;zr)] ' 

where Dy = BuBji — (BB T )ij is a positive-definite matrix. In a weak sense (when one is only 
interested in statistics of the process), one can speak of an equivalence between SDEs and Fokker- 
Planck equations. As we shall see below, this correspondence is the cornerstone of the proposed 
numerical approach: the pdf can be obtained by simulating the motion of stochastic particles, i.e. 
Eq. (fTU|) . In other words, real particles are replaced by stochastic particles which, if the model 
is suitable, reproduce the same statistics as the real ones. Indeed, in many problems of practical 
concern, the dimension of the reduced state vector is large and properties of the coefficients A 
and B make the direct solution of the above PDE (Partial Differential Equation), i.e. the Fokkcr 
Planck equation, numerically difficult. Instead, it is more appropriate to calculate p r (t; z r ) (or any 



G 



moment extracted from it) from Eq. (I10[) . Practically, this is done by resorting to a dynamical 
Monte Carlo method, that is by simulating a large number N of independent realisations of Z r (t). 
Then, at each time step, the discrete pdf, p r N (t; z r ), can be computed from the set of N independent 
samples {Z r ' n (t)} as 

1 N 

f N (t- z) = -]T 5(z{ - Z[< n (t)) 5(4 - Z r 2 > n (t)) . . . 5(z r d - Z r d n (t)), (12) 

n— 1 

where d is the dimension of the reduced state vector and n stands for the sample index. The 
question to be answered is: in what sense does the ensemble {Z r ' n (t)}, from which the discrete pdf 
p r N (t;z r ) is extracted, represent the underlying pdf p r (t;z r )7 The answer to this question can be 
given in a weak sense, that is when convergence in distribution is ensured, i.e. p r ^(t', z r ) — > p r (t; z r ) 
when N — > oo. 

A sequence of random variables {A™} converges in distribution to A if and only if, for any 
bounded continuous function g on E, one has (g(X n )) — > (g(X)) when n — > oo. In calculations, 
the mathematical expectations ( • ) is estimated by the ensemble average ( • ) jv over TV independent 
samples. The law of large numbers tells us that (g(X)) jy is an unbiased estimation of (g(X)}, that 
is (g(X))x — > (g(X)) when N — ¥ oo. Then, according to the central limit theorem, the error 
e N — (.9(A)) at — (<?P0)> which is a random variable, converges in distribution to a Gaussian 
random variable of zero mean and standard deviation a[g(X)]/yN provided that the variance of 
g{X), a 2 [g{X)], is finite. 

In the following, the PDF approach shall always be understood as the numerical solution of the 
set of SDEs equivalent, in a weak sense as explained above, to the corresponding Fokker-Planck 
equation. 

2.2 Stochastic integrals and calculus 

Stochastic differential equations require a strict mathematical treatment. The mathematical speci- 
ficities of SDEs have far reaching consequences for the derivation of accurate numerical schemes, 
cf. Section [5. 31 Some basic explanations are now given to highlight the important points that are 
needed for practical purposes. As a matter of fact, Eq. ([TO]) is just a shorthand notation for 

Zl(t) = Zl(t )+ [ Ai(s,Z r (s))ds+ f B l3 (s,Z r (s))dW 3 ( S ), (13) 

Jt Jt a 

where the first integral on the RHS (Right-Hand Side) is a classical Riemann-Stieltjes one. In 
the second integral, integration is performed with a measure, dSST(t), that has non-conventional 
properties. The Wiener process can be defined |14] as the only stochastic process with independent 
Gaussian increments and continuous trajectories (an increment, over a time step dt, is defined as 
dWj(t) — Wj(t + dt) — Wj(t)). The Wiener process has the following properties [Ml IT5]: 

(i) the trajectories of Wj(t) are continuous, yet nowhere differentiable (even on small time 
intervals, Wj(t) fluctuates enormously), 

(ii) each increment is a Gaussian random variable: (dWj(t) 2p+1 ) = for the odd moments, 
(dWj(t) 2 ) = dt and (dWj(t) 2p ) = o(dt), V p > 1, for the even moments. Increments over 
small time steps are stationary and independent, (dWj (t)) = 0, V t, and (dWj (t) dWj (t 1 )) = 
with t^t', 

(iii) the trajectories are of unbounded variation in every finite interval. 

The last property is the reason why the treatment of stochastic integrals differs from that of 
classical (Riemann-Stieltjes) ones (the Wiener process is not of finite variation [15]). Without 
going too deep into mathematical details, property (iii) simply implies that speaking of a stochastic 
integral without specifying in what sense it is considered lacks rigour (in this work, all stochastic 
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integrals will be considered in the Ito sense). In classical integration, the limit of the following 

sum (<r fc G [tk,tk+i]) 



,t N 

/ B ij (s,Z r (s))dW j (s)= lim Y, B ^' Lr i T k))( W ^k + x)-W j {t k )), (14) 

J to N ^ +O °t / 

should be independent of the choice of T k . This is not true in the above integral, because of 
property (hi) 1 . As a consequence, a choice has to be made for the sake of consistency. Two 
main choices (there exist others) are encountered in the literature, the Ito and the Stratonovich 
definitions. In the Ito definition, r k — t k and the following limit is under consideration 

N 

lim VB^fe.Z^tlK^fe+il-^fe)). (15) 

N— >+oo * — * 
k=l 

This choice has a major drawback, i.e. the rules of ordinary differential calculus are no longer 
valid. However, this drawback is balanced by the zero mean and isometry properties which are of 
great help when deriving weak numerical schemes, see Section [5.3.11 

ti 

X{s)dW{s)) = 0, 

( 16 ) 

X(s)dW(s) / Y(s)dW(s))= (X(s)Y(s)} ds. 

t Jt! Jt! 

where ( • ) is the mathematical expectation (to < t\ < ^ < ^3, X(t) and Y(t) are two stochastic 
processes). These properties no longer hold in the case of the Stratonovich interpretation but 
the rules of ordinary differential calculus remain valid. In the Stratonovich interpretation of 
the stochastic integral, the basic idea is to choose r k as the middle point of the intervals, i.e. 
2Tfc = tk + tje+i- There are, as a matter of fact, several possible choices, the most commonly 
encountered in the mathematical literature being 

t 

B ii («,Z r (a))odW^(s) = 



N 



lim ^-[By^Z'Xtfc)) +B ij (t k+1 ,Z r (t k+1 ))](W j (t k+1 ) - W s (t k )), 

N— >+oo £ — » Z 



(17) 



where o indicates that the stochastic integral is treated in the Stratonovich sense. 

The distinction between the Ito and the Stratonovich interpretations is critical, especially when 
ensuring consistency in the derivation of weak numerical schemes. There is actually an equivalence 
between the two interpretations. In can be shown [Ml [16] that Eq. ([TO]) written in the Stratonovich 
sense 

dZT(t) = Ai{t,Z r (t))dt + B^^Z^t^odWjit), (18) 
is equivalent to the following SDE, written in the Ito sense, 

dzm = Mt,z r (t))dt+ B kJ (t,z r (t)) dBij{ !; zr{t)) dt 

oz k (19) 
+ B ij (t,Z r (t))dW j (t). 

This result can explain why, in some works, computations performed with identical models can 
lead to contradictory results (the difference between the results is a drift term which is, most of 
the time, not negligible). Let us stress, once again, that even though the purpose of the present 



1 Hore, the different modes of convergence of the random variable -the stochastic integral- are not dealt with. 
For further information, see 1141 1161 for example. 
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paper is not to present mathematical subtleties, a good understanding of stochastic calculus is 
needed when deriving weak numerical schemes for SDEs encountered in fluid mechanics. 

As in the present work the Ito interpretation is retained, let us briefly present the basics of 
stochastic calculus. It has been mentioned that, in the frame of the Ito interpretation, the rules of 
ordinary differential calculus are no longer valid. As a matter a fact, this non-trivial consequence 
can be understood by going back to the properties of the Wiener process. For the second order 
moment of the increments of W(t), one has (dW(t) 2 ) = dt. This non trivial result (in classical 
calculus one would expect a second order term) is inherent to the nature of the Wiener process (it 
is a non-diffcrentiable process). As a consequence, the rules of classical calculus must be modified 
when considering terms of at least order 2. This is the well-known Ito formula. For any stochastic 
process X(f) which verifies the following SDE 

dXi(t) = Ai(t,X(t))^ + Sy(t,X(t))d^(i), (20) 

the SDE verified by any smooth function /(i, X(t)) is 

df(t,x(t)) = ^L(t,x(t))dt+dXi(t)^-(t,x(t)) 

at oxi , . 

i a 2 f (21) 

+ -(BB T )jj(t, X(i))-^— ^— (i, ~K(t))dt, 

where the last term on the RHS is a new term with respect to classical differential calculus. 

The main tools of the general formalism and the basics of stochastic calculus have now been 
introduced. This short presentation is a brief summary and a reader willing to derive weak 
numerical schemes for SDEs may refer to Refs. [14l [16] for the mathematical background, Ref. [9] 
for the physical background and more importantly Ref. |17j for the derivation of numerical schemes. 
A key point, which is not developed here (the main point being that stochastic integrals require 
a careful treatment), is the derivation of stochastic Taylor series, a tool which is needed when 
attempting to develop weak numerical schemes for SDEs. There exists a comprehensive book on 
these techniques, see Ref. |17) . 



3 PDF models 

The purpose of this section is to put forward the system of equations which is solved in the present 
mean- field/PDF approach, and to reveal the existing link between the physics of the problem and 
the tools that were presented in Section [2J Once this is done, attention shall be focused on a 
central specific theoretical issue: the treatment of the multiscale character of the dynamics of the 
discrete particles, cf. Section @] 

3.1 Derivation of a PDF model 

If the trajectories of a pair of particles (a fluid particle and a discrete particle) can be modelled 
by writing a system of equations given by Eq. (| 10[) , two issues must be solved: 

(i) what is the dimension of the reduced state vector (what variables should be retained)? 

(ii) and what is the form of the coefficients (the drift vector and the diffusion matrix)? 

A physical answer can be given to issue (i) when there is a separation of (time) scales. Let dt 
be the reference timescale at which the physical phenomena are observed. The separation of 
scales is defined in terms of slow and fast variables. A slow variable is a variable whose integral 
timescale, T, is much larger than dt and vice-versa for a fast variable whose integral timescale is 
r, i.e. r -C dt -C T '. The answer comes from the application of ideas known in synergetics, the 
so-called slaving principle 18]. In this equilibrium hypothesis, the fast variables are assumed to 
relax 'very rapidly' to their equilibrium values which can be expressed as a function of the values 
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taken by the slow modes. A practical application of this principle is the fast- variable elimination 
technique where the fast variables are replaced by models which represent their equilibrium values 
and usually involve white- noise terms [5]. The fast- variable elimination technique can be used in 
the derivation of one-point PDF models for single-phase turbulent flows. 

In this work, answers to issues (i) and (ii) are given separately for the fluid and the particles in 
the form of one-point PDF models. In the one-point PDF approach for the fluid (see Section l3.1.ip . 
single-phase turbulence is under consideration whereas in one-point pdf models for the discrete 
particles, the fluid-particle mixture is under investigation with known properties for the fluid (see 
Section I5.1.2[) . The two-point description briefly introduced in Section [5] will be addressed in 
Section OS 

3.1.1 One-point PDF models for single-phase turbulent flows 

When the Reynolds number is sufficiently large, for a reference time scale dt in the inertial range 
(t„ <C dt <C Tl, where Tl is the integral Lagrangian timescale and t v is the Kolmogorov timescale), 
the Kolmogorov theory [J] tells us that, for Lagrangian statistics, the covariance matrix of velocity 
increments has the form 

(dU u (t) dU f j (t + dt)) = C (e) dt % , (22) 

where e(t, x) is the local instantaneous energy dissipation rate and Co is a constant. Equation (|22l) 
implies that one has for the autocorrelation coefficients Rjj = 1 — [(Co d£)/(2T£)\ ~ 1 (velocity) 
and Ra = t^/Tl <C 1 (acceleration) 8 . This shows that, for dt in the inertial range, the velocity 
of a fluid particle, U/(t), is a slow variable whereas the acceleration, Af(t), is a fast variable. 
This suggests, according to the slaving principle, that Af(t) should be eliminated and replaced by 
a function of the slow modes, position and velocity, i.e. Z r (i) = {x*(i), Uy(t)}. The Kolmogorov 
theory gives answers to issues (i) and (ii) : the dimension of the reduced state vector is 2 and the 
diffusion matrix should be given by Hy = W Cq (e) §ij . The use of a SDE is not enforced by Eq. 
(I22p but the linear dependence in time of the covariance matrix is a strong indication. For further 
justifications concerning the use of SDEs for modelling purposes, see Ref. [8]. 

Using different arguments, several researchers |13( I19j have shown that a Langevin equation 
model for single-phase turbulence is 

1 Q/p) 

dU fti (t) = ^ dt + GijiUfj - {Uj))dt + ^C (e)dWi(t), (23) 

where P(t, x) is the local instantaneous pressure field. All mean fields, i.e. (P), (U) and (e) 
are evaluated at time t for x = x^-(t). The return-to-equilibrium matrix, G^, depends on mean 
fields and is usually written Gy = — 5y/Tx, + G°- where is a timescale given by l/T^ = 
(1/2 + 3Co/4)(e)/A; (k is the turbulent kinetic energy). The anisotropy matrix, G°-, also depends 
on mean fields only and can take different forms |13j . 
Equation (|23[) has two noteworthy properties: 

- the coefficients of this SDE depend not only on time and Z r (t) (as in Eq. $10\i ) but also on 
the expected values of functionals of the state vector. This dependence of the coefficients 
has important consequences, not only for the mathematical formalism, but for the numer- 
ical algorithm, see Section 15.11 These equations are often called Mac-Kean SDEs in the 
mathematical literature. 

- the model is not self-contained since external fields are needed to compute the drift vector and 
the diffusion matrix, i.e. the state vector should rather be written Z r (t) — {x/(i), U/(t), s(t)}. 
In complete models, a specific SDE is written for e(t) = e(t,x.f(i)) and the mean pressure 
field is given by a Poisson equation, see Ref. [T3] for example. 

3.1.2 One-point PDF models for the discrete particles 

Let us suppose that the discrete particles are moving in a turbulent flow whose mean fields are 
known (only one-way coupling is considered for the moment, i.e. the presence of particles does not 
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modulate the turbulence). These mean fields are most commonly the two first velocity moments, 
(U(t,x)) and (U(i,x) ® U(i,x)), and (P(i,x)) and (e(f,x)). 

In the case of discrete particles, the choice of a suitable state vector is more difficult than in 
the fluid case since there are no general results indicating a clear separation of scales. However, an 
extension of Kolmogorov theory [8] 120) shows that a linear dependence for the covariance matrix 
of the increments of the fluid velocity seen, i.e. (dU 8) i(t) dU s j(t + dt)), can be obtained under 
some specific hypotheses, for dt in the inertial range. Without being a formal proof, this result 
suggests to include the fluid velocity seen in the state vector, Z r (t) = {x p (t), XJ p (t), U s (t)}, and 
to model the increments of U s (i) with a Langevin equation. This choice differs from the one 
inherent to kinetic models where Z r (t) — {x p (t), U p (t)} [HI [21]. The existing correspondence 
between these two approaches has been discussed elsewhere |S]: it can be shown that including 
the fluid velocity seen in the state vector presents several advantages from the modelling point of 
view. Furthermore, kinetic models can be retrieved from the Langevin models for U s (£) [5]. Issue 
(i) has now been addressed, the dimension of the state vector is 3. Let us move to issue (ii), that 
is to write a SDE for the increments of U s (t). 

From a physical point of view, the problem of modelling particle dispersion (i.e. deriving 
a model for JJ s (t), see Eq. (JTJ) is more complicated than the diffusion one (fluid particles, cf. 
Section [3.1.1j) since two additional physical mechanisms have to be accounted for: particle inertia 
characterised by the timescale t p and external force fields (gravity in our case), Eq. (JTJ. Two 
main approaches can be found in the literature: 

- approaches based on paths (trajectories). A two-step construction is considered: a La- 
grangian step and an Eulerian step. The Lagrangian step corresponds to the trajectory, 
over a time interval dt, of a fluid particle located at time t in the vicinity of the discrete 
particle (this step is directly given by Eq. (|23|)). The Eulerian step corresponds to a spatial 
correction which gives, from the location of the fluid particle at t + dt, the fluid velocity 
seen by the discrete particle at time t + dt. This modelling point of view has two major 
drawbacks: it leads to an artificial decrease of the integral timescale of U s (i) (denoted T£ i 
in the present paper) and there is no clear separation between the effects of t p and g [8] . 

- approaches based on the physical effects [22]. A two-step construction is also considered by 
decoupling the two physical mechanisms: the first step corresponds to the effects of t p in the 
absence of external forces (in that case T£ i varies between two limit values, Tg -the integral 
Eulerian time scale- when r p — > +oo and Tl when t p — > 0). The second step corresponds 
to the effects of gravity alone which induce a mean drift and result in a decorrelation of 
U s (i) with respect to U/(i). This effect is called the crossing trajectory effect (CTE) and 
is related to the mean relative velocity (U r ) = (U p — U s ). 

In the present work, the derivation of a model for XJ S (t) is carried out by resorting to an approach 
based on the physical effects where the influence of the first step is neglected. This assumption 
allows us to extend Kolmogorov theory since the increments of the fluid velocity seen are only 
governed by mean quantities [5]. 

Assuming, for the sake of simplicity, that the mean drift (the mean relative velocity (U r )) 
is aligned with one of the coordinate axes (the general case is discussed in Ref. [8]), it can be 
shown [20] that a possible model for the increments of the fluid velocity seen is (the summation 
convention by repeated indices does not apply to the third and fourth term on the RHS) 



dU s ,i(t) = -J^ Ld t+ i(Upj) ~ (Uj)) ^ dt 



l~ (U s ,t - (Ui)) dt + J (e) (cobik/k + ^(kk/k - 1)J dW t {t). 



(24) 



The CTE has been modelled by changing the timescale, compared to the fluid case, in the drift 
term (third term on the RHS) and by adding a mean drift term (second term on the RHS). The 
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time scale is modified according to Csanady's analysis [53] 



n„ = T L {l + P ^) V \ (25) 

where /3j = /3, if axis 1 is aligned with the mean drift, with j3 — Tl/Te, and in the transversal 
directions (axes labelled 2 and 3) ft = 2/3. In the diffusion matrix, a new kinetic energy has been 
introduced 

where u(t,x) = U(i,x) - (U) and 6, = T L /T£ -. 
Equation (|24[) has two noteworthy properties: 



- it is consistent, by construction, with Eq. (|23[) when r p — > 0, that is when the discrete 
particles behave like fluid particles, 

- it is a Mac-Kean SDE even though the mean fields of the fluid are known (they are given 
by solving RANS equations). Indeed, it is necessary to compute the mean velocity of the 
particles (U p ) to calculate not only the mean drift term (second term on the RHS) but also 
the integral time scale of U s (i), T L ((U s ) is also needed for the computation of this time 
scale). 

Moreover, it must be emphasised that Eq. (|24|) is a possible choice among others and that 
the exact form of a Langevin equation for U s (i) still remains an open issue, see for exemple 
Refs. [24j [25] for models suited for homogeneous turbulence. There exists an alternative to Eq. 
([24} in the literature [26] where the coefficients are slightly different (the drift vector and the 
diffusion matrix) , the main difference being the form of the mean drift term which is written in 
terms of instantaneous velocities rather that mean velocities, i.e. (U p j — U s j)(d(Ui) /dxj). This 
difference has been discussed elsewhere [§]. This form of the mean drift term does not change 
the methodology which is presented in the rest of the paper, but it modifies the structure of the 
system of SDEs, i.e. U s (t) depends explicitly on the particle velocity, U p (t). This matter will be 
discussed in Section [B] 

When two-way coupling occurs, that is when the mass of particles per unit volume of fluid 
is sufficient to influence the characteristics of turbulence, Eq. (|24[) can be supplemented by an 
acceleration term A p _,. s which accounts for the influence of the discrete particles on the statistics 
of the fluid velocity sampled along the trajectory of a discrete particle, i.e. 

dU„,i(t) = A s>i dt + A p ^ Sil dt + B Stij dWj(t), (27) 

where the drift vector A s and the diffusion matrix B s are directly given by Eq. (|24p . 

For Ap_> s , the underlying force corresponds to the exchange of momentum between the fluid 
and the particles (drag force). The acceleration acting on the fluid element surrounding a discrete 
particle can be obtained as the sum of all elementary accelerations (due to the neighbouring 
particles) [8] 

A^ s = - Qpft)Us ~ Up , (28) 
a f p f t p 

i.e., at the discrete particle location x p , the elementary acceleration (U p — U s )/r p is multiplied 
by the probable mass of particles divided by the probable mass of fluid (since the total force is 
distributed only on the fluid phase). In Eq. (|28p . it is implicitly assumed that all particles under 
consideration have the same acceleration. Moreover, a. fit, x) and a p (t, x) represent the probability 
of presence of the fluid and the particles, respectively (a/ + a p = 1). 
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The complete set of SDEs which describes the one-point dynamical behaviour of the discrete 
particles is 

dx p s(t) = U p ,i dt, 

dUp i ( t ) = U °* ~ V*i dt + g . dti (29) 
t p 

dU Sil {t) = A Sii dt + A p ^ Sil dt + B s<ij dWj(t), 

where A s and B s are calculated by resorting to Eqs. (|24 |) - ([26| . This set of SDEs is under 
investigation in the present paper and the numerical methods needed to solve it will be discussed 
in Section [5] 

3.2 Mean-field/PDF approach 

As specified in the Introduction, a mean- field/PDF approach is adopted here for the computations 
of polydispersed turbulent two-phase flows, i.e., the fluid is described by solving RANS equations 
whereas the dynamics of the discrete particles are simulated with SDEs. The general formalism 
presented in Section [2] and the models presented in Sections 13.1.11 and 13.1.21 are sufficient to derive 
the mean-field/PDF model which is used in the present work. As a matter of fact, the set of 
SDEs, from which the pdf of the discrete phase can be extracted, has already been presented, 
cf. Eqs. (|29|) . Only the derivation of the set of mean-field (RANS) equations for the continuous 
phase has to be put forward. This derivation can be done using different techniques 27, 28] [29]. 
An interesting technique, which is in line with the tools introduced in Section [2j is to resort to 
a two-point description. This new approach hardly provides further physical information but it 
allows both the fluid and the particles to be described under the same formalism. This description 
is briefly outlined here, supplementary information is given in Appendix [XJ 

3.2.1 PDF models for polydispersed turbulent two-phase flows 

The path which is adopted is to gather the preceding results that have just been presented for the 
time increments of the fluid velocity seen along discrete particle trajectories, Eq. (|2"9"1) . and for the 
time increments of the fluid velocity along fluid particle trajectories, Eq. (1231) . The system of SDEs 
is, however, supplemented by one term (an acceleration) Aj,_^ which reflects the influence of the 
discrete particles on the fluid. The time rate of change of Z r (i) = {x/(t), U/(t), x p (i), TJ p (t), U s (t)} 
is given by 

' dxf.i(t) = Uf t i dt, 

dUf t i(t) = Af ti dt + Ap-+f ti dt + Bf ti j dWj(t), 
< dx Pti (t) = U Pti dt, (30) 

dU Pji (t) = ((Z7 Sj j - U p j)/t p ) dt + g { dt, 
k dU Stl (t) = A 3ti dt + A p ^ Sil dt + B„ tij dW-(t), 

where the drift vector Af and the diffusion matrix B f are given by Eq. ([23]). The form of A p _j./ is 
discussed in Appendix [A] By assuming that the trajectories of a pair of particles can be obtained 
in such a way, i.e. Eqs. p0[) . one should be aware that several assumptions have been made. 
Further explanations on these assumptions are given in Refs. [S] and |29) . 

3.2.2 Mean- field/PDF model 

It can be shown that the mean- field (RANS) equations for the fluid can be extracted from Eqs. ([30]) . 
see Appendix [A] The complete set of equations is given in Table [TJ Depending on the closures 
chosen for the return-to-equilibrium matrix, Gy, different RSM equations can be obtained. In 
practical computations, these mean-field equations are supplemented with a PDE for (e) obtained 
by the same path (as mentioned in Section ["3.1.1[ a SDE can be written for e(t)). 

The mean- field/PDF model has now been presented, that is the mean field equations for 
the fluid and the set of SDEs for the discrete particles. In practical computations, the mean 
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fluid properties are computed with classical finite volume techniques whereas the dynamics of the 
stochastic particles are calculated by resorting to stochastic calculus, i.e. by integrating in time 
a set of SDEs (weak numerical schemes). The finite volume algorithms are well known and no 
discussion on these methods is given here. However, some specific issues related to the computation 
of the source terms in Eqs. (jl 29[) and (|130j) will be addressed in Section 16.11 These source terms 
are computed from the stochastic particles (the numerical solution of the SDEs) and, as we shall 
see in Section 15.11 information has to be exchanged between the stochastic particles and the mesh 
on which the mean-fluid properties are calculated. Before carrying on with the derivation of weak 
numerical schemes for the set of SDEs, a very important issue has to be put forward: the multiscale 
character of the model for the discrete phase. 

The discussions to come on the multiscale character of the set of SDEs, Section [H and on the 
derivation of weak numerical schemes for these equations, Section [5. 31 are presented in the case of 
one-way coupling. This is not due to a limitation of the method but it simply reflects the current 
status of the present work. Indeed, the methodology, which is presented in Sections 0] and [5j 
remains valid in the case of two-way coupling. The presentation of the general algorithm, Section 
15 . 1 1 and the discussion on issues related to projection and averaging in particle-mesh method, 
Section 15. 2\ are put forward in the case of two-way coupling. Specific issues related to two-way 
coupling will be discussed in Section [6j namely the computations of the source terms in the set 
of PDEs (RANS equations) and the extension of the material presented in Sections 0] and 15.31 to 
two-way coupling. 

4 Multiscale properties of the SDEs 

There are three different timescales describing the dynamics of the discrete particles, cf. Eq. (|2T))) : 
dt, the timescale at which the process is observed, T£ i the integral timescale of the fluid velocity 
seen, XJ s (t), and t p the particle relaxation time. Once again, it must be recalled that these SDEs 
have a physical meaning only in the case where dt <C T£ i and dt <C r p . What happens if one of 
these conditions or both are not verified? It is in fact possible to show that the system of SDEs 
converges towards several limit cases which are consistent with the physics. The mathematical 
details are not given here, see Ref. [30 for further explanations. Nevertheless, the fast-variable 
elimination technique is now presented in a simple case in order to help the reader understand the 
form of two of the limit cases of system (f29|) . 

4.1 Fast-variable elimination technique 

Let us consider a model problem to illustrate this technique. A good historical example is the 
treatment of Brownian motion. There exist two points of view to address this problem: Einstein's 
point of view where only position is retained in the state vector, Z(i) = {x(i)}, and Langevin's 
approach where the state vector is composed by position and velocity, Z(t) = {x(i), U(i)}. Let us 
consider one-dimensional Brownian motion for the sake of simplicity. Langevin's model reads 

dx(t)=U(t)dt, 

dU(t) = -(U(t)/T)dt + BdW(t), [ ' 

this system being valid when dt <C T (T is the integral time scale of U(t), B is the diffusion 
coefficient). What happens if this condition is not verified, i.e. when T — > 0? The velocity, U(t), 
becomes a fast variable and, according to the slaving principle, it should be eliminated and its 
influence should be expressed as a function of the slow modes (here position, x(t)). This new 
description should then be consistent with Einstein's model, dx(t) = \J2D dW(t) where D is the 
diffusion coefficient. 

By applying the rules of stochastic calculus (see Section [2T2T) . one can show that 

i r* 

U(t)~BTr)(t) with r/{t) = -exp(-</T) / exp(s/T) dW(s), (32) 
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where the influence of initial conditions have been neglected (one integrates from — oo). Here, rj(t) 
is a Gaussian random variable [15] (it is a stochastic integral of a deterministic function) with 
(r)(t)) = (zero mean property) and (isometry property) 

(V(t)v(t')) = -!= exp(-|t - t'\/T) — - > 5(t - t>). (33) 

21 T-> 

Therefore, rj(t) is a Gaussian white noise and consequently dW(t) = 77(f) dt. The limit system for 
Langevin's model, Eqs. pip , when T — ¥ and BT remains finite, is then given by 

dx(t) = V2D dW(t), (34) 

where D — (BT) 2 /2. The velocity, whose integral timescale becomes zero and whose variance 
becomes infinite (Gaussian white noise) has been eliminated, that is the process is observed at 
a timescale dt large in comparison to the velocity fluctuations (characteristic time scale, T). Its 
influence is, however, left in the diffusion coefficient of the reduced model (Einstein's model) 
provided that the product BT remains finite. As we shall see shortly in Section |4~2| this model 
problem is going to be helpful to understand two of the four limit cases to come. 

The above analysis of the model problem has been presented in the continuous sense, that is 
when both variables are continuous functions of time. In practical applications, i.e. in numerical 
computations or in experiments, both variables would have been observed at discrete times. In 
numerical calculations, one computes the solutions at discrete times where the time step corre- 
sponds to the observation timescale dt. In experiments, the variables of interest are measured at a 
given sampling frequency / and one has dt = 1//. The development of the numerics is slightly an- 
ticipated and the limit case T — ¥ is now investigated in the discrete sense, that is when T <C At. 
A good understanding of the subtle difference between the continuous case and the discrete case 
will be helpful in the development of weak numerical schemes to come, see Section [5.3. II 

In the discrete sense, the limit case does not correspond to T — > but rather to T < At. 
Applying the rules of stochastic calculus, the solution to system (j3"Tj) at t = to + At is given by 

f x(t) = x(t ) + U(to) T[l- exp(-At/T)] + I x {t), 
\U{t) = U(t )exp(-At/T) + I u (t), 

where the stochastic integrals are defined as 



I x {t)=BT [ dW{s) - BTexp(-t/T) [ exp(s/T) dW(s) 
Iu(t) = Bexp(-t/T) I exp(s/T)dW(s). 

J tg 



(36) 



Here, (I x ,Iu) is a vector composed by two dependent, centred (zero mean property, see Section 
2.2[) Gaussian random variables (Gaussian since I x and Ijj are stochastic integrals of deterministic 
functions 15 ). It can be shown that a centred Gaussian vector can be expressed as the product of 
two matrices by resorting to the Choleski algorithm: these matrices are the covariance matrix and 
a vector composed of independent standard Gaussian random variables (A/"(0, 1), i.e. zero mean 
and variance equal to unity). This decomposition is well suited to numerical applications since a 
set of independent standard Gaussian random variables can easily be generated on computers by 
using suitable random number generators. By applying the Choleski algorithm to (I x ,Iu), one 
can write 

i x (t) = ((i x iu)/J{i^) zu + J(i*)-(i x iu)V(ih)t x , 

v / (37) 



15 



where £ x and £j/ are two independent standard Gaussian random variables. The components of 
the covariance matrix are given by 

{I 2 ) = {B T) 2 {At - T[l - cxp(- At/T)} [3 - exp(- M/T)]/2}, 
{I 2 ) = B 2 T[1 - exp(-2At/T)]/2, (38) 
(I x I v ) = {BT[1 - exp(-At/T)]} 2 /2. 

Therefore, in the discrete sense, the limit system to Langevin's model becomes 

Pf 



x{t) = x{t Q ) + U{t )T + BT 




(39) 



/ B2T 



In the discrete sense, the velocity U{t) does not 'disappear' (in the continuous case, it becomes 
Gaussian white noise). This result is physically sound since the velocity is only observed at time 
steps which are large compared to its memory (integral timescale). Finally, the above results are 
consistent with the observation of Einstein, that is, for long diffusion times, one has {x{t) 2 ) ~ 
{BT) 2 t. 

4.2 Limit cases 

From now on, the summation rule by repeating indices is dropped to avoid confusion, as in Eq. 
(|2"41 for example. The system of SDEs describing the dynamics of the discrete particles reads 
{from now on B a ^j is denoted Bij for the sake of simplicity) 



' dx P a(t) = U v<i dt, 
dU p At) = — {U s a - U p ,i) dt + A dt, 

dU sA {t) = - — U Sti dt + C t dt + ^2 B i:i dW 3 {t), 



(40) 



where Ci is a term that includes all mean contributions: the mean pressure gradient, —(d{P) jdxi)j pf, 
the mean drift term, {{U p j) — (Uj))(d(Ui) jdxj), and the mean part of the return-to-equilibrium 
term, ([/,) /T£ As explained above, two-way coupling is left out of the present analysis. A4 is an 
acceleration (gravity in the present work, but it can be extended for practical reasons to the case 
of other external force fields). Once again, system (|40p has a physical meaning only in the case 
where dt <C T£ i and dt <C r v . When these conditions are not satisfied, it is possible to show that, 
in the continuous sense (time and all coefficients are continuous functions which can go to zero), 
the system converges towards several limit systems |30j . 
Case 1: when r p — > 0, the particles behave as fluid particles and one has 

dx p j{t) = U Pl i dt, 

syM.o.u(40) > ) U pM = U s ,i(t), (4l) 



T„->0 



dU s 4t) = - — ^dt-^{U s ^{U l ))dt + ^/Q^)dW l {t), 
Pf axi ±l 



that is, the model is consistent with a known turbulent fluid PDF model [11] as explained in 
Section f3. 1.21 This shows that the model is a coherent generalisation of the fluid one, which can 
be recovered as a limit case. 

Case 2: when T£ . — > and BijT£ ■ —> est, the fluid velocity seen becomes a fast variable. It is 
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then eliminated and one can write 



™» {BijTl ; ->cst) 

system (M : > < 



dU p ,i(t) =—({Ui) - Up,,) dt + A dt+ 

T v (42) 



E 



13 L < 1 dWj(t). 



This result can be understood from the model problem of the previous Section. U s (i) has been 
eliminated but its influence is left in the diffusion coefficient BijT£ i /r p . In this case, the equations 
are equivalent to a Fokker-Planck model for particles of significant inertia. 

Case 3: When t p , i — > and at the same time B^T^ i — > est , the fluid velocity seen becomes 
a fast variable and the discrete particles behave as fluid particles. It can be shown that 

{dxpjt) = (U l )dt + A l dt+ 
^(B^dWM (43) 

We retrieve a pure diffusive behaviour, that is the equations of Brownian motion, cf. Section [4.11 
Case 4: at last, when T£ i — > with no condition on B^T^ i; the velocity of the fluid seen is no 
longer random and the system becomes deterministic. The flow is laminar and it can be proven 
that 

dx Pii (t) = U Pli dt, 

system gDJ — ► { dU Pti (t) = —{{Ui) - U Pji ) dt + A, dt, (44) 

U.,i(t) = (Ui). 

Limit cases 1 to 3 reflect the multiscale character of the problem. When the timescales go to 
zero (with a condition on their products with the coefficients of the diffusion matrix) , a hierarchy 
of stochastic differential systems is obtained. Moreover, the elimination of the fast variables (the 
velocities U p (t) and U s (t)) does not mean that these variables do not (physically) exist anymore: 
they simply become Gaussian white noise. As we shall see in Section [5.3. 11 both ~U p (t) and U s (i) 
become independent Gaussian random variables, in the discrete sense, in limit case (iii) (in limit 
case (ii), only U s (i) becomes a Gaussian random variable). These results are in line with the 
previous model problem, cf. Eq. (|39| in Section |4~T1 

The existence of limit systems is a key point in the development of weak numerical schemes 
to integrate in time the set of SDEs describing the dynamics of the discrete particles, i.e. Eqs. 
(|40p. As we shall see in the next Section, in numerical computations, dt the observation timescale 
of the process, is the time step. A suitable weak numerical scheme should therefore be consistent 
with all limit cases since, as we shall see, it is not possible to control the time step to enforce 
the conditions necessary for the validity of Eqs. ([40l) . Before we carry on to the time integration 
of Eqs. (|40p . let us give a general overview of the numerical procedure which is needed to solve 
the whole set of equations (the mean-field (RANS) equations for the fluid and the SDEs for the 
discrete particles). 



5 Numerical approach 

The mean- field/PDF model used in the present paper for practical computations has now been 
given, see Table [TJ It consists in a set of PDEs describing the dynamics of mean-fluid quantities 
and a set of SDEs from which the joint pdf of the variables of interest for the discrete particles can 
be extracted. In this approach, the numerical solution is obtained by resorting to a hybrid method 
where the mean- fluid properties are computed by solving the mean- field (RANS) equations with a 
classical finite volume procedure whereas the local instantaneous properties of the discrete particles 
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are determined by solving the set of SDEs, Eqs. ([2"9"]l . Therefore, the mean fluid properties are 
computed on a mesh whereas the statistics of the discrete phase are calculated from particles 
moving in the computational domain. 

A closer look at the equation system presented in Table [T] shows that the set of equations has 
the following properties: 

(i) in the set of PDEs from which the mean-fluid properties are computed, mean fields involving 
discrete particles properties are needed in order to compute the source terms, see Eqs. (|129|) 
and 11301) . 

(ii) in the set of SDEs, the knowledge of statistical moments (for the fluid and the discrete 
particles), such as mean values and variances, at the locations of the stochastic particles, 
is required in order to compute the time evolution of the discrete particle properties and, 
thus, of the statistics derived from them. Indeed, from a mathematical standpoint, the set 
of SDEs can formally be written as (from now on, the notation is slightly changed: Z(t) 
denotes the state vector for the discrete particles, i.e. the superscript r is dropped for the 
sake of clarity) 

dZ[t) = A(t,Z(t),p(t;z),Y(t))dt + (r(t,Z(t),p(t;z),Y(t))dW(t), (45) 

where p(t;z) stands for the pdf of Z(t) and Y(t) represents external mean fields, i.e. the 
fluid mean fields defined at particle locations. For each time t, p(t; z) has to be calculated 
in order to compute the coefficients of the SDE: in the present approach, the pdf (or any 
necessary moment extracted from it) is computed out of all stochastic particles that are 
tracked (all values taken by Z(i)), cf. Section I2TT1 Thus, as mentioned in Sections 13.1. II and 
13.1.21 a kind of integro-differential equation, no longer local in the space of Z(i), is obtained; 
it is called a Mac-Kean SDE and is inherently difficult to solve [31] [32]. In other words, we 
are dealing with a system where particles interact weakly (or indirectly) through the mean 
fields that they create (it is the leading idea of one-point PDF models). 

As far as property (ii) is concerned, in practice, probabilistic expectations of particle properties 
at a given point are approximated by spatial averages over nearby particles, i.e. the statistics 
extracted from the stochastic particles (which arc needed to compute the coefficients of system 
([2"9"))) are not calculated for each particle (this would cost too much CPU time) but are evaluated 
at each cell centre of the mesh (generated for the solution of the set of PDEs) following a given 
spatial average (averaging operator). These moments can then be evaluated for each particle by 
interpolation or projection (projection operator). The same projection operator is used to compute 
the statistics of the fluid at the locations of the stochastic particles. The source terms in the PDEs, 
cf. property (i), are directly computed by resorting to the averaged particle quantities. 

This is the principle of particle-mesh methods: exchange data between particles and mesh 
points. In the present work, the main advantage of such methods is of course the reduction of 
CPU time but the use of projection and averaging operators has some drawbacks, i.e. it creates 
additional numerical errors and, in the general algorithm, each particle has to be located in the 
mesh. 

5.1 Particle-mesh methods 

Historically, particle-mesh methods have been widely used in other areas of physics like the dy- 
namics of plasmas, astrophysical simulations, electrostatics, etc. [33]. In these applications, the 
system of equations which is solved is deterministic and the mesh is uniform, most of the time for 
unbounded domains or bounded domains with periodic boundary conditions. These models are 
often referred to as the particle-in-cell (PIC) approach [531 [511, 155] . 

In fluid dynamics, apart from calculations of dispersed two-phase flows, particle-mesh meth- 
ods are mainly encountered in computations of single-phase turbulence with stand-alone particle 
models [36, 37, 38 , in hybrid particle/field models [39], in calculations of turbulent reactive flows 
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with PDF methods [TU] , and in numerical simulations of non- Newtonian polymeric fluids [21] ■ In 
these applications, the situation is different from the classical PIC approach. Indeed, the systems 
of equations are stochastic and the domains are bounded with various boundary conditions. Let us 
present the main lines of the particle-mesh algorithm used in the present work. It is implemented 
in the ESTET software [40]. 

5.1.1 General algorithm 

Let {Yl x l} stand for the set of the fluid mean fields at the different mesh points and let {Y^* 1 } 
be the fluid mean fields interpolated at particle locations. Let {Z^'} denote the set of variables 
attached to the stochastic particles and {Z^l} the set of statistics, defined at cell centres, extracted 
from {Z^- 1 }. Time is discretised with a constant time step At = t n+ i—t n and space with a uniform 
mesh of cell size Ax. 

The first step (operator F) of the algorithm is to solve the PDEs describing the fluid, 

{YM}(t„) and {ZM}(f n ) A {Y^}(t n+1 ). (46) 

The F operator corresponds to a classical finite volume RANS solver and it gives the evolution in 
time of the statistical moments of the fluid (the particle properties are needed to compute source 
terms when two-way coupling is accounted for). 

The second step (projection, operator P) consists in calculating mean-particle properties and 
mean-fluid properties at particle locations, 

{zH}(tn) and {YW}( tn ) A {Z^}(t n ) and {Y^}^), (47) 
Then, the stochastic differential system can be integrated in time (operator T), 

{zW}(f„) and {YW}(g A {zW}(t n+1 ). (48) 

Finally, from the new computed set of variables, at particle locations, new statistical moments 
are evaluated at cell centres, 

{zW}(t n+1 ) A {ZW}(t n+1 ), (49) 

and so on. The general algorithm is therefore defined by iterating the four operators, F — > P — > 
T — > A. The purpose of the present section, Section [5] is to discuss different implementation 
aspects of the general algorithm and more especially issues related to: 

(i) the specificity of the treatment of the averaging and projection operators, Section I5T121 

(ii) the time integration of the SDEs (operator T), i.e. the determination of a suitable weak 
numerical scheme for system (|29l) which is consistent with the multiscale character of the 
physics (asymptotic cases), Section 1531 

Before discussing these issues, some last clarifications are given on the nature of the numerical 
errors generated by the particle-mesh algorithm sketched above. 

5.1.2 Numerical errors in particle-mesh methods 

The numerical particle-mesh solution of evolution equations like Eq. (1451) involves several kinds of 
errors. These errors have been described in the context of PDF methods for turbulent reactive 
flows [37] [41]. The overall error of the PDF computation (P — s- T — >• A) can be separated into a 
deterministic and a random part, the former involving the bias, spatial and temporal discretisation 
errors. Numerical errors occur due to: 

(i) spatial discretisation, represented by a typical mesh size Ax, 

(ii) finite temporal resolution, determined by the time step At, 
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(iii) the use of a finite number of particles, both in the whole domain (N) and per cell (N pc ); 
these are further decomposed into the statistical error (of zero average) and the bias. 

The spatial discretisation error (i) is akin to the classical error in numerical methods for solving 
PDEs and depends on the mesh size. In the present approach (Mac-Kean SDEs), it is inherent 
to the use of projection and averaging operators (A and P). This numerical error does not 
occur when the coefficients of the SDEs depend only on local values of the state vector. This is 
the case, for instance, in numerical computations of stand-alone PDF methods for homogeneous 
turbulence [37]. To date, only Xu & Pope [41] have addressed this issue for non-linear SDEs and 
have found some characteristics of the error for an infinite number of particles per cell. Note that 
the spatial discretisation error occurs also in grid- free methods, such as SPH [42], and is, in this 
case, directly related to the smoothing parameter (kernel size). 

The temporal discretisation error (ii) is basically the same as in any numerical method for 
solving the time evolution of the solutions to deterministic ODEs or PDEs. Numerical schemes 
for the SDEs (operator T) must be developed and analysed with care due to the specificity of 
stochastic calculus, cf. Section I5~3I 

The statistical error, which is inherent to any Monte Carlo method, is due to the use of a finite 
number of particles per cell (samples) to compute the statistics and is proportional to the inverse 
of the square root of N, according to the central limit theorem. In specific applications (e.g., [36]). 
the coefficient of proportionality can be reduced considerably when appropriate variance reduction 
techniques (VRT) are applied [37JII3]. 

The bias is the difference between the mean value of a quantity for a finite number (N) of 
particles and the mean value for infinitely many particles (all other parameters being unchanged), 
i.e. for any random variable Z 

1 N 

B Z (N) = ((Z) N )-(Z) with {Z)n = ^Y. Z ^ ( 5 °) 

i=l 

where Z l stand for the different sample values of Z . The bias is thus a deterministic error, 
important for non-linear stochastic models [37] HT] . The issue is worth explaining with a simple 
example. Consider a random variable X with a certain law (pdf), say standard Gaussian, i.e. 
(X) = and (X 2 ) = 1. We note X € W(0, 1). For the mean value (X) N computed out of N 
samples, the central limit theorem gives (for N sufficiently large) 



(X} N = {X) oa + C£/VN, (51) 

where £ £ A/"(0, 1) and C is a proportionality constant for the statistical error; obviously (X)^ = 
{X} = 0. Consider next a function Y = (X) 2 . Now Yoq = 0, but for a finite number N of 
samples Y N = C 2 £ 2 /N and after averaging (Y N ) = C 2 /N. The bias B Y (N) = (Y N ) - Y x is thus 
proportional to N^ 1 . In more general terms [37], for Y = g((X)) we have = g((X)oo), and 
the development into a Taylor series, accounting for Eq. (|5ip . yields 

Y N = g({X) N )=g({X) O0 + -%=£ 



N ' / 

Cn' C 2 n" { a'" x 



/JV 2N s V^ 3/ 
After averaging the above, the bias is computed as 

*rW=C 2 ^ + o(j^) . (53) 

It depends on the local second derivative of g and is proportional to iV _1 . In a general case of 
random fields, the bias interplays with the spatial error because of the kernel estimation which is 
applied to compute averages [13l [44] . 
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5.2 Averaging and projection operators 

We recall that, in the numerical solution process of Eq. (|45p . moments of Z(i), like (Z) and (ZiZj) 
need to be extracted from the particle data. A correct computation of these quantities is crucial for 
the overall numerical solution, since the moments are put back into the SDE and serve to advance 
the particle properties to the next time step of the simulation, see Section [5. 1.11 These ingredients 
of the algorithm, i.e. the computation of mean fields (or averaging) and their interpolation at (or 
projection to) particle locations are well known in the PIC approach [3311311133]. For those particle 
models (deterministic equations on a regular mesh), optimum averaging and projection schemes 
have been worked out. In the present case, a new problem is addressed: stochastic models with 
boundary conditions typical of fluid mechanics. Consequently, some new important numerical 
features appear. 

(i) The first specificity is related to the computation of statistics: in most applications, we need 
not only the mean values, but also (at least) second-order moments present in the evolution 
equations. These moments are usually position-dependent (non- homogeneous in space). 

(ii) A second specificity here (often present in applications in fluid mechanics) is that the com- 
putational domain is bounded and the associated mesh is non-uniform; as argued further, 
adequate boundary conditions may affect the computation of statistics. 

Here, attention is focused specifically on the errors due to the exchange of information between 
particles and the mesh, that is how mean fields (usually the first and second-order moments) are 
computed and projected at particle locations, the main issue being to investigate whether classical 
techniques already used in particle simulations are also suitable for our present particle-mesh 
problem. 



5.2.1 Averaging operators, weighting functions 

In order to introduce the numerical issues related to the exchange of information between the 
particles and the mesh, let us first discuss the difference between the ensemble mean (expected 
value) and the spatial average. For the case of a deterministic function $(x), the spatial average 
( ■ )a (with a characteristic smoothing length A) in the cell centred at x™, i = 1, . . . , J, can be 
thought of as the integral: 

($ [l] ) A =J $(x) «J(x - xW) dbt, (54) 

where $W = $(xW) and w is a given weighting function (smoothing kernel) satisfying J w)(x) <ix = 
1. For a random field x) with a pdf /$(£, x; ^P), the mean at corresponds to the probabilistic 
expectation, i.e. 

/+oo 
tt/*(t,xM;*)tftt, (55) 
-oo 

where \t is a sample-space variable of $. In classical Monte Carlo methods, the pdf /$(t,xM; ty) 
at each point x^ is approximated by using a number N of independent samples of the random 
variable, say & n >,n = 1, ...,N. In other words, a set of variables is attached to every 
particle, located at x < - n - ) in the computational domain (NB: the superscript convention helps to 
distinguish between particles ^ and cell- related quantities W ) . 

The mathematical expectation, in Eq. (1551) . is computed exactly at x^ , whereas a in 

Eq. (|54|) represents the spatial average centred on x^ ; both are equal in the spatially-homogeneous 
case only. However, space-dependent moments cannot be calculated exactly from a Monte Carlo 
estimation, given a finite number of particles used in the simulation. In practice, under a local ho- 
mogeneity assumption, expectations ($W) are approximated by (local) spatial averages ($>^)n,a, 
based on a discrete particle set: 

N 

~ jv,a = $( "M x(n) - * W ) • (56) 

71=1 
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This expression, also known as kernel estimate |13j . is derived from Eq. (|54p as a quadrature 
formula; the weights w are the discrete equivalents of w, X^=i w ( K< " n ^ ~ x '*') = 1- 

A generalisation of the above discussion to centred moments of any order is straightforward. 
Here, the explicit formulation is given for the variance because of its importance in further con- 
siderations, cf. Sections 15.2.21 to 15.2.41 By analogy to Eq. (fS"5")) . the local value (at x^) of the 
variance cr| of the random field <£>(£, x) is 

/+oo 
(*-<$W)) 2 /»(*;xW,t)d*. (57) 
-oo 

As an extension of formula (|56l) for the mean value, the expression for the variance <j\ with the 
use of spatial averaging becomes 

N 

(oi)W * (4)S,a = E( $W - (*»^(x (n) - xW), (58) 

n=l 

where ($) stands either for a cell average or is interpolated at xW. From the algorithmic 

standpoint, Eq. (|5"5)) has some disadvantages. Firstly, the double pass over particles (the mean 
has to be computed in advance) results in some computational overhead. Secondly and more 
important, for first order weighting functions (cf. Sections 15.2.21 and I5.2.3[) . it implies a risky 
extrapolation of mean fields at the locations outside the computational domain. We propose to get 
around the difficulty by computing the central moments directly from the standard (non-centred) 
moments; for example, the variance of a random variable Q satisfies ((Q — (Q)) 2 ) = (Q 2 ) — (Q) 2 - 
Yet, some precaution is needed since such an expression for the variance can not be guaranteed 
to remain always non-negative because of round-off errors. 

As far as the choice of the weighting function is concerned, two methods are widely used 
and are under investigation in the following: the NGP (Nearest-Grid-Point) and CIC (Cloud- 
in-Cell) methods. They correspond to weighting functions of zero (constant) and first (linear) 
order, respectively. These methods have been thoroughly discussed [33] . as mentioned above, for 
deterministic systems solved on uniform meshes and in the particular case of unbounded domains or 
bounded domains with periodic boundary conditions. In the present work, we address the problem 
for the numerical solutions of specific stochastic systems (Mac-Kean SDEs) on non-uniform meshes 
for bounded domains and non-periodic boundary conditions. 

In the NGP method, particle properties are associated with the centre of the cell containing 
the particle (for a uniform mesh, it is the grid point nearest to the particle, hence the name), and 
the weighting function is top-hat (or piecewise constant, Fig. QJ,), that is ui(x^™^ — x^) = 1/iVL. 
for particle n in cell i and otherwise (Np C is the number of particles in cell i). The NGP average 
is thus found from the sum over all N^ c particles in a given cell i 

1 Kc 

< $w > = ^rE $(n) - ( 59 ) 
J V „=i 

In the CIC method, for a uniform mesh in D spatial dimensions, the weighting function is piecewise- 
linear (Fig. QJ>): 

D 1 

»(x) = [|r«{l-|xj|/Ai,0}, (60) 

3 = 1 3 

where Aj is the width of the cell in direction j. The particle is thus regarded not as a single point 
but rather as a linear distribution of properties: a cloud centred at xW, with a width of 2Aj. 
The CIC method is less local than NGP: the average at a given cell centre is computed not only 
from the particles located in this very cell, but also from those in the neighbouring cells. 

Higher order, but less local, weighting functions can be used, for instance a piecewise-quadratic 
polynomial (Fig. [TJ:) or a cubic spline as in Smoothed Particle Hydrodynamics (SPH) [45]; quartic 
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and quintic splines as well as Gaussian kernels are also quite popular |421 146] . An alternative to 
kernel estimators, cf. Eq. (|56p are somewhat more costly (especially in 2D/3D) formulae with 
least-squares or local least-squares approximations, as well as cross- validated splines [TU1 144) . 

The projection of averaged (cell) values onto the particle locations is an interpolation procedure, 
in a sense akin to averaging, Eq. (|56|) . with W and M replaced by each other (it is the reverse 
operation). The consistency of averaging and projection steps has been a serious concern in PIC 
applications. Indeed, it has been demonstrated (see [33], §5.2.4) that, in the case of a system of 
charged particles moving in an electric field generated by themselves, unphysical forces may appear 
if the averaging scheme is not of order equal to (acceptably also higher than) the projection scheme. 
Those authors have also stated that the CIC method performs better than NGP. 

5.2.2 CIC statistics on a non-uniform mesh 

Although the usual presentation of CIC formulae is made for a uniform mesh, in typical appli- 
cations of the present mean-field/PDF approach for polydispersed turbulent two-phase flows, a 
non- uniform mesh may be of advantage (wall-bounded flows). A number of difficulties arise in 
this generalisation, depending whether the mean density computed from the particle masses or 
the mean of a variable attached to the particles is sought. For the sake of simplicity, we will illus- 
trate the issue in a ID setting; extension to 2D and 3D is straightforward through the Cartesian 
products. 

First, consider the computation of fluid density. A particle located at x^ n > where x^ < x^ < 
x^ l+1 \ cf. Fig. Ufa), will contribute to the mean density at xW and x^ +1 ^. Here, we treat the 
particle as a "cloud" (or a linear distribution of mass) centred at x^ n ' and stretched on the 
interval Sxi — min(Axi, Axj+i). With this assumption, the particle mass is contained between 
Xq 1 ' and Xq +1 ' and does not contribute to other cell averages (three possible locations of a 
"stretched" particle are shown in Fig. [2^). The fraction of the cloud that belongs to the cell 
[i + 1], R x = (x^ + Sxi/2 — Xq )/Sxi, adds to the average at x^ l+1 \ whereas (1 — R x ) adds to the 
cell average at x". For a boundary cell, the whole mass of particles located close to the boundary 
(between the boundary and the centre of the cell closest to it) is attributed to the centre of the 
cell next to the boundary. Fig. [3] presents 2D computation results of the r.m.s. density (r.m.s. 
deviation from the mean particle density) a 10*10 mesh, both uniform and non- uniform, using the 
standard NGP averaging and the CIC method described above. For the uniform mesh, particle 
locations in the domain are generated randomly from the uniform random distribution. The 
non-uniform mesh has been generated using random numbers from a uniform distribution on the 
interval (0.5Aa;, 1.5Aa;); in this case, particle locations have been generated deterministically with 
a constant number density. To reduce the statistical error of the r.m.s. density deviation, it has 
been computed as an average of four different runs with different seeding values for pseudorandom 
number generator. As expected, the statistical error is higher for the NGP average and varies as 

— 1/2 

N pc . We note that although the r.m.s. density on the non-uniform mesh is higher, there is no 
systematic error for the CIC computation which confirms the correctness of the method above. 

Next, consider the computation of cell-averaged values for any variable $ assigned to the 
particles (a typical example is the calculation of the mean velocity at a given point from the set 
of particle instantaneous velocities). In this case, both the particle masses and the values of $ 
attached to the particles come into play (it has been mentioned in Appendix [S] that, in the frame 
of the present formalism, one has to resort to mass- weighted expected values). The first idea 
could be to use the method described above for the calculation of the mean density. However, in 
a basic test of a linear deterministic function and particles distributed uniformly on the interval, 
this method is readily shown not to retrieve given values of particle variables at cell centres for 
a non-uniform grid. Fig. [2th) illustrates this point where cell i is half the size of its neighbours. 
To compute the average value of $ at the particles labelled 2, 3 and 4 contribute with the 

respective weights of 1, 3/4, and 1/4, 

2 = 1($ + A$) + ^($ - A$) + ~($ - 3A$). (61) 
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For this simple example, it is readily checked that depends on A$ (it should not) and 

that ^ $ - The reason for this behaviour can be traced into the formula of particle weights 

which is not symmetric with respect to x^ +1 ^. Indeed, the same behaviour can be noticed when 
non-symmetric CIC weight functions are used though they seem to be a natural generalisation 
of Eq. (|60|) . This point is also illustrated in Fig. [2^b) that includes one such basis function w(x) 
decreasing linearly from a maximum at the considered cell centre, here a;[ l+1 l, to zero at the 
neighbouring cell centres, here x^ and x^ l+2 ^. Using these weighting functions, we readily find 
again 

I (* [ml > = ^(*o + A3) + ^($ - A*) + i(* - 3A$), (62) 

that is (<£>[ t+1 ]) =/: $ 0l a result that would explicitly depend upon the gradient of the function 
and the size of the mesh. Instead, we propose the following method that works correctly for 
particle variables on a non- uniform mesh. We generalise Eq. ([50)1 by taking A as min{Aa;i_i, Ai J 
to preserve the symmetry of w{x). This expression gives the correct result for the example of 
Fig. [2fb). An alternative to the above is the use of a two-stage algorithm [44], satisfactory but 
arguably more time-consuming. Again, a 2D computation has been performed for the r.m.s. of a 
variable assigned to particles from a deterministic linear profile in space with the same procedure as 
above for particle density computations. Both the standard NGP averaging and the CIC method 
described above have been used. Results shown in Fig. 0] confirm the advantages of CIC averaging 
in this case. To continue, let us now remove the deterministic assumption for $ and consider 
averaging of a random variable assigned to particles 

$(«) = $(a;(»)) = n»(iW) + s(a;(™))£ , (63) 

where £ <E Af(0, 1), m(x) = do + a\X is a linear function and s{x) = A (where A denotes the 
average mesh size) to study the effect of the spurious variance resulting from the NGP averaging. 
A numerical test has been performed again with the same methodology as described above. Fig. [5] 
shows the normalised r.m.s. of this linear random function (i.e. the square root of its variance 
computed at the cell centres divided by the prescribed r.m.s.) computed out of particle locations. 
A systematic error is readily noticed for NGP statistics, unless aiA <C 0{1). 

5.2.3 CIC statistics with boundary conditions 

In wall-bounded flow applications, the CIC method has to be modified so that suitable boundary 
conditions (BC) are properly accounted for. Dreeben and Pope [44] describe their application of 
CIC statistics in the PDF method but the treatment of flow boundaries is not reported there. In 
Fig. |51 a schematic plot shows clearly why the CIC averaging requires some care in presence of flow 
boundaries (here walls). The NGP statistics are local (in a cell) and they can thus be computed 
in border cells without any change. However, for the CIC method this is not the case. As it 
transpires from Fig. [6] the CIC average of the particle density in a border cell, computed with 
the weight wi(x), incorrectly gives a smaller value than the density in the neighbouring internal 
cell, computed with W2{x). The obvious reason is that less particles contribute to the mass in the 
border cell. The problem differs depending upon whether we consider the particle density field or 
the mean of any variable $ (such as velocity) assigned to particles. To compute the mean density, 
for the particles located between the cell centre and the boundary, it is sufficient to assign all 
their mass to the centre of this boundary cell. For other particle variables, this treatment is not 
applicable: even for a linear function <S>(x), the higher spatial accuracy (compared to NGP) of 
CIC is lost. 

A working remedy to this situation is the addition of "ghost" or "mirror" particles outside the 
actual computational domain in order to compensate for the incorrect CIC computation at the 
cell centres close to the boundaries. The idea of these "ghost" particles is known in the context of 
the SPH simulation [46]. The very presence of ghost particles with masses equal to their "hosts" 
allows for a correct CIC density computation in boundary cells. Next, to compute CIC statistics 
of any variable $ attached to the particles, values $' of the variables are to be determined also for 
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the mirror ones, as they enter Eq. (1561) with a corresponding weight. The procedure is relatively 
straightforward for the CIC averaging of given functions (either deterministic or random) where the 
value of the function at a mirror particle location is known. However, the main interest for using 
ghost particles is in actual particle simulations where precisely these functions are unknown. The 
values of the variables attached to the mirror particles are to be determined directly from those of 
the host particles through the application of relevant BCs. For example, if $ stands for a particle 
velocity component, then the no-slip, impermeable wall implies $' = — <f> (Fig. [6]) so that at the 
boundary ($) = 0. Generalisation is possible for more complex, yet still block-structured, 2D and 
3D geometries. Other types of possible boundary conditions in the PDF method for turbulent flow 
computations, e.g. where particle boundary conditions are determined from a physical reasoning 
for the near- wall region |36j . are also readily implemented this way. 



5.2.4 Test case: space-dependent SDE 



Let us now analyse the behaviour of a simple but nevertheless realistic example of the generic 
SDE, Eq. (|4"5|) . where emphasis is put on the averaging and projection operators. Both the NGP 
and CIC techniques are going to be used with the suitable modifications presented in Sections 
15.2.21 and 15.2.31 We consider a one-dimensional Ornstein-Uhlenbeck process on < x < 1 interval 



d$(i) = - 



$(i) - m(x) 
f 



dt 



2s 2 (x) 



T 



dW(t), 



(64) 



with a fixed timescale T; the space-dependence of the mean value m(x) and the variance s 2 (x) is 
essential to study the spatial discretisation error, cf. Section [5.1.2l In Eq. (|64|) . there is no physical 
coupling in x; it would occur if convection were added (e.g. dx — /($) dt) or if non-local operators 
were present (e.g. aV 2 ($)dt terms on the RHS). Here, only a "numerical" coupling exists; it is 
due to the approximate computation of ensemble averages through spatial averages, Eq. (I56[) . 

For a given time step, At = £;+].— ti, trajectories of the stochastic process can be integrated 
analytically [51 117). thus avoiding any temporal discretisation error; this yields (the location x is 
a parameter) 



= exp(-^f) + m(x) 



l-exp(-^f) 



! s 2 (x) 



, 2 At 
1 - exp( — — ) 



(65) 



where £ is a standard Gaussian random variable, £ 6 A/"(0, 1). 

The actual test case consists in using the computed values of ($) and a\ in Eq. (|65|) in place 
of m(x) and s 2 (x), respectively. The computed profiles (at different times) for a quadratic mean 
m(x) and a constant variance s 2 {x) are presented in Fig. [71 Qualitatively, a typical behaviour of 
CIC averaging is to modify the shape of the mean ($>)(t,x) (ultimately towards a linear profile) 
and to produce an increase of the variance in the centre of the computational interval. On the 
contrary, the NGP mean value remains basically unchanged, while the variances in separate cells 
become increasingly uncorrelated in time. Indeed, in the CIC method, neighbouring cells are 
linked through the averaging procedure, thus the profile of the variance stays relatively smooth. 
The NGP computation of Eq. ([M]) is local and as a result, separate cells become independent from 
one another, cf. the upper right plot in Fig. [7J These numerical outcomes raise the question of the 
existence of a bias (cf. Section l5.1.2l) . A detailed numerical study of this issue has been performed 
by analysing the temporal behaviour of the r.m.s. of <f> averaged over all cells, denoted by u$(t) 
and given by 



<7*(t) = y5>M(f), 
i=l 



(66) 



where erM (t) are the r.m.s. values obtained in cell i. Indeed, for the NGP method, calculations done 
in each cell are uncorrelated and ct^ (t), i = 1, . . . , I, can be regarded as independent samples. The 
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simplest case to be discussed now is that of constant mean m(x) and variance s 2 (x) in Eq. (|64p. 
Fig. [SJa) shows how c$(f) normalised by its initial value, evolves in time. The procedure, in 
the NGP case, is sensitive to the number of particles per cell. In subsequent time steps, the 
recomputed variance decreases. The decrease is faster for a smaller number of particles per cell. 
This phenomenon is identified to be related to the bias and is due to the finite value of N pc . 
Fig. [8]Jb) clearly indicates that there is a bias; however, computations of the variance of c$(f) 
(results presented as error bars in the same figure) show that the statistical fluctuations of this 
quantity are significant. With this precaution in mind, the bias (resulting from several different 
runs for several number of cells, etc.) has been plotted as a function of the number of particles 
per cell N pc . Results in Fig. Hta) indicate that the slope for different runs is indeed close to the 
theoretical prediction of the bias, -B$ ~ (plotted as a dotted line). The computed probability 
distributions of 0$ are shown in Fig.^b). 

5.2.5 NGP or CIC in practical computations? 

In Fig. the centred second-order moment, i.e. the variance, has been computed using the non- 
centred moments in order to avoid problems with the correct statement of boundary conditions 
which are necessary for CIC averages. It is readily seen that the computed profiles of ($)(x) 
and <J 2 {x) tend to become linear (both are fixed to the respective m(x) and s{x) values at the 
boundaries x — and x — 1). In Fig. HI it is seen that the NGP computation generates a non-zero 
spurious variance. 

The comparison between the NGP and CIC methods for stochastic processes reveals that none 
of the results are entirely satisfactory. Compared to classical deterministic particle systems, where 
the CIC has a known advantage over NGP, new features have been revealed. NGP is a local 
and robust method which does not require specific developments. It respects given mean profiles 
but leads to spurious variances. The spurious variances are shown to be related to a bias which 
decreases linearly with the number of particles. On the contrary, when using CIC, developments 
are needed to account for non-uniform meshes and for boundary conditions. It is a less local 
(higher order) method than NGP but it is more complicated to implement in the general case. 
The results obtained in a prototypical SDE, cf. Section lS.2.41 show that the CIC method does not 
present advantages in accuracy. Furthermore, in the particular case considered (no convection), 
the CIC technique performs globally worse than the NGP one. Indeed, CIC does not preserve the 
mean and shows a systematic error in the variance. 

As a consequence, all practical computations performed in the rest of the present paper (cf. 
Section [7|) will be based upon the NGP technique. 

5.3 Accurate schemes for SDE integration 

It has now been explained how information is exchanged between the discrete particles and the 
mesh (operators P and A, cf. Section lip.l.ll) . We recall that this is done here in the case of one-way 
coupling and the extension of the numerical schemes, presented in this Subsection, to the case of 
two-way coupling will be discussed in Section [6] 

As shown in Sections 13.1.21 and 13.21 and as explained at the beginning this Section, cf. Eq. 
(|45p . the SDEs reproducing the dynamics of the discrete particles are Mac-Kean SDEs since the 
coefficients depend not only on the state vector but also on expected values of functions of Z(t). 
In particle-mesh methods, as explained above, quantities such as (/(Z)) (where / is a linear or 
quadratic function of Z(t) in our problem) are extracted from the particle data and evaluated at 
grid points. In Section 15.21 the difficulties inherent to the projection and averaging procedures 
have been detailed. 

Attention is now focused on the time integration of the set of SDEs (operator T). The de- 
velopment of a suitable weak numerical scheme for the time integration of SDEs is a much more 
difficult task than the corresponding one for ODEs. Indeed, SDEs do not obey the rules of classi- 
cal differential calculus, see Section l2~2l and one has to rely on the theory of stochastic processes 
[15] . In that sense, particular attention has to be payed to the problem of consistency between 
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discretised equations and the original continuous set of SDEs. It is recalled that, in the present 
paper, Ito's calculus is adopted and therefore all SDEs are written in the ltd' sense, see Section 
12.21 This choice has no physical motivation: Ito's calculus is very convenient in the development 
of weak numerical schemes for SDEs because of the zero mean and isometry properties, cf. Eqs 
CU). 

An essential preliminary is to clearly frame the development of the weak numerical scheme 
in the general methodology that has been followed here. We propose, indeed, to describe the 
problem of turbulent polydispersed two-phase flows within an engineering context but with a 
rigorous treatment of the multiscale character which is a distinctive feature of these flows, cf. 
Section [4] The model, Eqs. (|40|). contains several characteristic timescales and this system of 
SDEs becomes stiff, from a mathematical point of view, whenever one of these timescales goes 
to zero. In those cases, various limit systems are obtained, see Section 21 which represent the 
asymptotic limits of the physical model. A proper treatment of the physics of the multiscale 
aspect imposes to put forward weak numerical schemes which are consistent with all asymptotic 
limits when the different timescales go to zero. 

It is worth emphasising that this corresponds to a practical concern. Indeed, in the numerical 
simulation of a complex flow, the timescales may be negligible (much smaller than the time step) 
in some areas of the computational domain. For example, in a wall-bounded flow, the integral 
timescale of the fluid velocity seen, T£ -, goes to zero when the distance to the wall decreases. 
Furthermore, when dealing with polydispersed particles, or with phenomena involving evaporation 
or combustion (when particle diameters decrease in time), one has often to handle a whole range 
of particle diameters (say from lfim to 100 — 200^m) and thus a whole range of values for r p (for 
example, 10~ 6 s < r p < 10~ 2 s). It would be inefficient to carry out computations with a time 
step limited by the smallest possible value of r p and/or T£ 4 . This is the very reason why in the 
simulations of particle dispersion in wall-bounded flows, based on discrete models, it is necessary 
to use different time steps for different classes of diameters (and thus of t p ) and to lower the time 
step in the wall region. 

As a consequence of the above discussion, the constraints which are required for a suitable 
weak numerical scheme, can now be summarised, considering both physical and numerical issues. 
Since a particle-mesh method is adopted here, the PDEs for the fluid are first solved and then 
the dynamics of the stochastic particles are computed (see Section [5.1.11) . thus, the scheme has 
to be explicit for the fluid mean fields. By choice, the scheme will also be explicit for the particle 
properties. Furthermore, the time step, which has to be the same for the integration of the PDEs 
and the SDEs, is imposed by stability conditions required by the finite volume algorithm solving 
the mean-field equations for the fluid. This implies that, since there is no possibility to control the 
time step when integrating the SDEs, the numerical scheme has to be unconditionally stable. At 
last, since particle localisation on a mesh (needed for projection and averaging, see Section [Oil is 
CPU-demanding, the numerical scheme should minimise these operations. 

The constraints required for a suitable weak numerical scheme are: 

(i) the numerical scheme must be explicit, stable and the number of calls to particle localisation 
(sub)routine has to be minimum, 

(ii) the numerical scheme must be consistent with all limit systems. 
5.3.1 Weak numerical schemes for SDEs 

In Section [21 the correspondence (in a weak sense) between a set of SDEs and a Fokker-Planck 
equation (for the associated law) has been established. In this work, weak numerical schemes 
shall be developed for Eqs. (j40]l . i.e. we are not interested in the exact trajectories of the process 
but instead in statistics (the pdf) extracted from the stochastic particles (the real particles are 
replaced by stochastic ones which should reproduce the same statistics). The numerical method 
proposed in this work is therefore nothing else than the simulation of an underlying pdf, or in other 
words, the equivalent Fokker-Planck equation is solved by simulating the trajectories of stochastic 
particles, that is by a dynamical Monte Carlo method. As briefly explained in Section [2~T| this 
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non-trivial numerical procedure, i.e. to resort to SDEs to solve a PDE, is well suited for PDEs 
with a large number of degrees of freedom. 

Since the Ito interpretation of stochastic integrals has been chosen, it is implicitly assumed, 
in the discretisation of the stochastic integral, that Bij should not anticipate the future, i.e. for 
each time step At — tk+i — tk, should be computed at t = t/.- As a result, classical numerical 
schemes for ODEs, for example Runge-Kutta schemes, can not be applied directly. More precisely, 
careless applications of such schemes for SDEs can introduce spurious drifts which may not be 
easy to detect. An illuminating example of this kind of error is illustrated in Ref. |47j . The key 
point here is that the numerical discretisation of the stochastic integral must be in line with its 
mathematical definition. 

Let Zi At (t) be a numerical approximation of Z(i) obtained with a uniform time discretisation, 
At. A numerical scheme of order m will converge, in a weak sense, if at time t — T (T is called 
the stopping time), for all sufficiently smooth functions /, there exists a constant C (function of 
T) such that 

su P |(/[Z(t)] -f[Z At (t)})\ < C(T)(At) m . (67) 

t<T 

Other convergence modes are possible, for example strong convergence in the mean-square sense, 
if one is interested in the exact trajectories of the process. It is fairly rare that this is the case for 
engineering problems. Indeed, in most engineering applications, one is mainly interested in the 
expected values (statistics) of functionals of the variables of interest. For further discussion on the 
convergence modes, see the book of Kloeden & Platen [T7] . 

5.3.2 Analytical solution to the system of SDEs 

In the present work, the weak numerical schemes, with the required features, are developed based 
on the analytical solution to Eqs. (|40|) with constant coefficients (independent of time), the main 
idea being to derive a numerical scheme by freezing the coefficients on the integration intervals. 
This methodology ensures stability and consistency with all limit systems: 

- stability because the form of the equations gives analytical solutions with exponentials of 
the type exp(— At/T) where T is one of the characteristic timescales (t p and T£ f ), 

- consistency with all limit systems by construction, since the schemes are based on an ana- 
lytical solution. 

Different techniques shall be used to derive first and second-order (in time) schemes from the 
analytical solutions with constant coefficients. A first-order scheme can be obtained by computing, 
at each time step, the variables on the basis of the analytical solutions (all coefficients are frozen 
at the beginning of the integration interval), i.e. a numerical scheme of the Euler kind is obtained. 
A second-order scheme can be derived by resorting to a predictor-corrector technique where the 
prediction step is the first-order scheme. 

Before presenting the weak numerical schemes, it is a prerequisite to give the analytical solu- 
tions to system (|40l) . with constant coefficients (in time). These solutions are obtained by resorting 
to Ito's calculus in combination with the method of the variation of the constant. For instance, 
for the fluid velocity seen, one seeks a solution of the form U S: i(t) — Hi(t) exp(— t/Tj), where Hi(t) 
is a stochastic process defined by {from now on the notation is slightly changed: T£ i is noted Ti 
for the sake of clarity in the complex formulae to come) 

dHi(t) = exp(t/Ti)[Ci dt + Bi dW^t)], (68) 

that is, by integration on a time interval [to, t] (At = t — to), 

U, ti {t) = U„,i{to) exp(-At/T 4 ) + C t Ti [1 - exp(-At/T)} 

ft (69) 
+B l exp(-t/T l ) / expis/TijdWiis), 

J to 
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where Bi = Bu since Bij is a diagonal matrix, cf. Eq. By proceeding in the same way 

for the other equations (position and velocity), the analytical solution is obtained for the entire 
system, cf. Table [5J The three stochastic integrals, Eqs. (| 134[) to (1136)) in Table [H are centred 
Gaussian processes (they are stochastic integrals of deterministic functions 15 ). These integrals 
are defined implicitly, but they can be simplified by integration by parts, cf. Table [2] As explained 
in Section 14. 1[ for the numerical representation of the stochastic integrals, the knowledge of the 
covariance matrix (second-order moments) is needed, see Table El Using the isometry property, 
see Section |2~2| the second-order moments, i.e. Eqs. (|140j) to ()145j) in Table EJ can be calculated. 
The analytical solutions are now known. Before presenting the first order scheme, let us verify 
that the analytical solution given by Tables [2] and [3] is consistent with the limit cases obtained in 
Section [OJ i.e. Eqs. (|4T1) to (144 ] , 

Limits systems of analytical solution: in limit case 1, where the discrete particles behave 
as fluid particles, the limit system is given by Eq. (|41l) . When t p — > 0, Eq. (|132|) becomes 

U p ,i(t) = U s<i (to) exp(-At/T,) + a T l exp(-At/T t ) + T z (t), (70) 

and for the stochastic integral I\(i), one has 

<I?(*)) -> ^[1 - cxp(-2Ai/T,)] = ( 7l 2 (t)), 

(Ti(t)7i(*)) ^(*)>- 

The last two equations indicate that Ti(t) — > Ji(t) when t p — > 0. By comparing Eq. ([70| to Eq. 
(fT33"l) with Ti(t) = 7i(t), the results of Eq. fl4l]) are retrieved, i.e. U p (t) = U s (t). 

In limit case 2, the fluid velocity seen, XJ s (t), is a fast variable which is eliminated. The results 
obtained in Table [5] and [3] with Ti — > and Bi Ti = est, give 

UpAt) = U p .i(t ) exp(-At/T p ) + [(Ui) + A t p ][1 - exp(-At/r p )] 



/ B 2 T 2 (72) 
+ W^ r i [l-cxp(-2Ai/r p )] g p , h 

where Q p ^ is a 7V(0, 1) vector (composed of independent standard Gaussian random variables) and 
we recall that (Ui) — (Ui(t,x p (t)). It can be rapidly verified, by applying Ito's calculus, that Eq. 
(1721 is the solution to system (|4"2")l when the coefficients are constant. 

In limit case 3, both the fluid velocity seen and the velocity of the discrete particles become 
rapid variables. When t p — > and Ti — > with Bi Ti = est, Eq. (|131|) becomes 

x P ,i{t) = x p ,i(t ) + [(Ui) + Ai t p ] At + ^B 2 T 2 At g s>i , (73) 

which is the solution to Eq. (|4"3"]) when the coefficients are constant (Q x> i is a Af(0, 1) vector). 

In limit case 4, when Ti — > (with no condition on B{Ti) the system becomes deterministic, 
the results are in agreement with Eq. (|gg|). When T t 0, Eqs. (TilTLj) to (TT^)) become 

r u sA (t) = 

U P}i {t) = U Pti (t Q ) exp(~At/T p ) + [(Ui) + A r p ][l - exp(-At/r p )], 
x Pti (t) = x P:l (t ) + t p [1 - exp(-At/Tp)]J7p,i(t ) 

+ [(Ui) + At r p ]{At - t p [1 - exp(-AVr p )]}, 

which is the analytical solution to system (|44|) when the coefficients are constant. 
5.3.3 Weak first order scheme 

The derivation of the weak first order scheme is now rather straightforward since the analytical 
solutions to system (|40|) with constant coefficients have been already calculated. Indeed, the Euler 
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scheme (which is a weak scheme of order 1 [IT]) is simply obtained by freezing the coefficients at 
the beginning of the time intervals At = [t n ,t n+ i]. Let Zf and Z™ +1 be the approximated values 
of Zi(t) at time t n and t n +i, respectively. The Euler scheme is then simply written by using the 
results of Tables [5] and [3] as shown in Table 2] Before showing that the scheme is consistent with 
all limit cases, some clarifications must be given. Here, the limit systems are considered in the 
discrete sense. The observation timescale dt has now become the time step At. The timescales t p 
and Ti do not go to zero, as in the continuous sense (Section |4]) , but their values, depending on 
the history of the particles, can be smaller or greater than At. The continuous limits, i.e. Eqs. 
([4"Tj) to (U3}, represent a mathematical limit, whereas in the discrete formulation, as we shall see 
just below, the limit systems correspond to a numerical solution where the ratios At/Ti and At/r p 
become large (the limit systems are obtained by Taylor expansions). 

In limit case 1, when t p — > in the continuous sense and r p <C At <C Ti in the discrete sense, 
the numerical scheme gives U^f 1 — U™f , see Table 01 which is consistent with the results of 
Section [O] 

In limit case 2, in the continuous sense Tj — > and B{Ti = est, that is the fluid velocity seen 
U s it) becomes a fast variable which is eliminated. In the discrete case, U s (t) is simply observed 
at a timescale which is great compared to its memory, that is Ti <C At <C r p , and the numerical 
scheme yields (see Table E} 



where (£/") = (£/i(t„, x™)). The fluid velocity seen becomes a Gaussian random variable, a result 
which is physically sound since XJ S (t) is observed at time steps which are greater than its memory. 
This result is in line with that of the model problem presented in Section 14.11 Furthermore, by 
Taylor expansion, it can be shown that the numerical scheme is consistent with Eq. (I72|) . 

In limit case 3, that is when 1 <C At/Ti and 1 <C At/r p (discrete case), one obtains for the 
velocity of the particles and for the fluid velocity seen (see Table 0J 



Once again, U p ^(t), and U s .i(t), which were eliminated in the continuous case, do not disappear. 
They become Gaussian random variables, a result which is physically sound since these two random 
variables are observed at time steps which are greater than their respective memories. Moreover, 
by Taylor expansion, one can show that the numerical scheme is consistent with Eq. (1731) . 

In limit case 4, Ti — 0, and the flow becomes laminar. It can be easily shown that the numerical 
scheme is consistent with Eqs. (f74"]) . For instance, one has for the fluid velocity U^f 1 = (J7™), cf. 



The previous results show that the Euler scheme presented in Table [4] is consistent with all 
limit cases. Therefore, the scheme gives numerical solutions which are physically sound, i.e. a 
consistent representation of the multiscale character of the model is obtained. 

5.3.4 Weak second order scheme 

Most of the time, dynamical Monte Carlo methods are used with first-order schemes only, for 
example in nuclear or particle physics. In those cases, the time-step value is not a very important 
factor, and attention is rather focused on obtaining accurate statistics. On the contrary, in in- 
dustrial fluid mechanics applications with complex geometries and strong inhomogeneities in the 
flow, a high-order accuracy in time can be critical in order to avoid prohibitively small time steps 
resulting in huge computational time. Such an example will be presented in Section 17.11 




(75) 




(76) 



Table H 
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From a formal point of view, weak high-order schemes for a set of SDEs can be derived for any 
desired accuracy, though this is much more complicated than for ODEs. Such high-order schemes 
are generally based on truncated stochastic Taylor expansions, see for example Refs. [TTJ H5] . 
These techniques can not be applied directly in our particular case since neither the unconditional 
stability nor the consistency in limit cases can be obtained. 

Property of the system of SDEs: the diffusion matrix of system (|29l) has a singular prop- 
erty of crucial importance here |30| . In the present case, this nine-dimensional matrix can be 
written, using block notation, as (we recall that Z(£) = (x p (t), TJ p (t), TJ s (t))) 



<r(t,Z(i)) 





B s (t )Xp (t)) 



(77) 



where each block represents a three-dimensional matrix. Indeed, from Eq. (I24|) . it can be no- 
ticed that B s i j depends only on time, position, x p , and the mean value of the relative velocity, 
(U r ). Therefore, the only variable of the state vector on which B s ^ depends is position, because 
(U r )(t, x p (t)) is a mean field. The fact that cry depends neither on U p nor on XJ S implies that 
quantities such that d<Jij/dzk are non-zero only when 1 < k < 3 and 6 < i, j < 9. For these values 
of k and j, one has Ukj = 0. 

Thus, the diffusion matrix cry has the following singular property 



\ ~* \ ~* d&ij 



0, Vi. 
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General idea: let us consider the following model problem 

dXi(t) =StipC(t))dt + 5^58«(X(t))dW^(i), (79) 

j 

where 03^ verifies property (|78p . It can be shown, for example by stochastic Taylor expansions 
[17) . that a predictor-corrector scheme of the type 

1 / 3 n 1 / x ( 8 °) 

x? + - (a? + ar 1 ) At + Y, 5 (®S + 

is a weak second-order scheme (2t^ +1 = %(X n+1 ), <B™ +1 = < By(X™+ 1 ), At = t n+1 - t„ and 
AWj = W" +1 — Wj l ). This result is true, once again, only when the diffusion matrix verifies 
property (f75|) . cf. Ref. [57] • If this property is not verified, the problem is more complex and 
other terms are needed to enforce second-order accuracy, see for example Talay [35]. Since the 
predictor step of the scheme above is the Euler scheme (already developed in Section [5. 3. 3[) . the 
remaining task consists in finding a suitable correction step which ensures the fulfilment of the 
constraints listed above. 

Derivation of the numerical scheme: how should the coefficients of the predictor step, 
and be computed? The main idea here is to generate a correction step based on 

the analytical solutions by considering that the acceleration terms vary linearly with time. This 
idea originates from considerations related to Taylor series expansions. The numerical solution 
obtained from the analytical solution with constant coefficients is an approximation of first-order 
accuracy. Mathematically, the solution is given in terms of the integral of acceleration terms. 
Thus, one can state that the solution based on the zero-th order (constant terms) development of 
the acceleration terms gives a first-order approximation in time. By analogy, it can be guessed that 
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approximating the acceleration terms by piecewise linear functions in time yields a second-order 
approximation in time. 

Let us introduce the following notation: U^f 1 and U^f 1 stand for the predicted velocities 
and If +1 and f p n+1 are the predicted time scales. The values of the fields related to the fluid 
taken at (t n +i, x™ +1 ) are denoted, for example, (f or (P n+1 ). As far as the computation of 
the mean fields extracted from the discrete particles are concerned, it is worth emphasising that 
none of them are computed at (t„+i, x™ +1 ), because the scheme would become implicit, i.e. fields 
such as the expected value of the particle velocity are computed from the predicted velocities. For 
example, one has 

(TT n + 1 \ 

a(t n+1 , x^ 1 ) = cr +1 = ^t+r + /(( u ; i+1 >. (u n+1 ). <^ l+1 »- (si) 

i 

Let us first consider the fluid velocity seen. The analytical solution to system when the 
coefficients are constant is, by applying the rules of Ito's calculus 

U s At) = U Sji {t ) exp(-At/T 2 ) + f Ci(s,x p )exp[(s - t)/T;] rfs+ 7i (t), (82) 

Jt a 

where the temporal coefficients (the timescales) are considered constant, while the term Ci(s,x p ) 
is retained in the integral. Following the previous ideas, let us suppose that Ci(s, x p ) varies linearly 
on the integration interval [to, t], that is (At = t — to) 

a-(s )Xp (s)) = C i (to,x p (t )) + ^[a(«,x p (t))-C i (to,x J) (to))](s-to). (83) 
By inserting Eq. (|53")) into Eq. integration gives 

U s>i (t) = U.,i(to) exp(-At/T 2 ) + [T i C i {t ^ p {t ))]A 2 (At,T i ) 

+ [T i C i (t,x p (t))]B 2 (At,T i )+ li {t), ( ' 

where the functions ^(At, x) and _B2(At,x) are given by (x is a positive real variable) 

( A 2 (At,x) = -exp(-At/x) + [1 - exp(-At/x)][At/x], 
1 B 2 (At,x) = 1 - [1 -exp(-At/x)][At/x\. 

Accounting for the time dependence of the coefficients, i.e. Tj, it is proposed to write the following 
correction step, which is in line with the treatment of the acceleration terms, 

Kt 1 =\ Ki exp(-At/I™) + \ U? fi exp(-At/Tf +1 ) 

+A 2 (At, Tf) [T?C?] + B 2 (At, fl l+1 ) [ff +1 C™ +1 ] +7™ +1 , 

where a consistent formulation for the stochastic integral 7™ +1 is needed. The same procedure is 
used, i.e. the diffusion matrix B,j is linearised and integration is carried out. The final expression 
is 

/ rpn+l 

7f +1 - y [Sf] 2 ^-[1 - exp(~2At/T;" +1 )] g M , (87) 

where Gi,i is the same 7V(0, 1) random variable used in the simulation of 7," in the Euler scheme 
and where B* is defined by 



l-exp(-2At/T™ +i ) 



/; 



A 2 {2 At, 77+ 1 ) J{B?y + B 2 (2 At, Tf +1 ) J (B? +1 V 
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Here, some explanations must be given. During integration, another step is necessary in order 
to achieve the closed form presented in Eq. (|88p . Indeed, two parts derive from the integration 
by parts carried out when Bij varies linearly. The first term is an analytical function, while the 
second term is still a stochastic integral, therefore the global integral can be written formally 
7™ +1 = 61+62- It has been considered that a projection of this second integral on the first 
remains of second-order accuracy for the global scheme. Therefore, the following hypothesis has 
been used, 8 2 ~ ((<5i(5 2 )/(<5 2 ))<5i. 

In the case of the velocity of the particles, the same approach followed for the fluid velocity 
seen is adopted. Let us start from the exact solution with constant coefficients for XJ p (t). By 
resorting to the rules of Ito's calculus, one can write 



U p ,i(t) = U p ,i(to)exp(-At/T p )+ 

1 f* 

— exp(-At/Tp) / exp(s/T p )[U s ,i(s) +TpAi(s,x p )]ds, 
t p J t 



(89) 



and by inserting Eq. (|82[) in the previous equation, one has 
U p ,i(t) = Up ti (t ) exp(-At/T p ) + U s ,i(t ) 9 t [cxp(-At/T l ) - exp(-At/r p )] + I\(t) 



H exp(-t/r p ) / exp(s/r p ) 



exp(-s/T l ) / Ci(u,x p ) exp(u/Ti)du + TpAi(s,x p ) 



ds. 



(90) 



Two deterministic integrals must be treated in Eq. (|90[) . A multiple one, involving Ci(t,x p ) and a 
simple one with the acceleration term Ai(t,x p ). Both integrals are handled as done previously for 
the fluid velocity seen, that is, it is assumed that both accelerations vary linearly on the integration 
interval, see for example Eq. (|83|) for Ci(t,x p ). By integration by parts of both integrals, one 
finds after some derivations 
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U p ,i(t) = Up yi (t )exp(-At/Tp) 

+ U„,i{to) 0<[exp(-Ai/2<) - exp(-At/r p )] 
+ [T l C l (t ,Xp(t ))]A 2c (Tp,T l ) + [T l C l (t,x p (t))]B 2c (Tp,T l ) 
+ [TpA t (t ,Xp(t ))}A 2 (At,Tp) + [TpA l (t,x p (t))]B 2 (At,Tp) 
+ Ti(t), 

where the functions C 2c (x,y), A 2c (x,y) and B 2c (x,y) are given by (x and y are two positive real 
variables) 

C 2c (x,y) = [y/(y-x)][exp(-At/y)-exp(-At/x)], 

A 2c {x, y) = - exp{-At/x) + [{x + y)/At] [1 - exp(-A*/a;)] 

-(l + y/At)C 2c (x,y), 

B 2c (x, y) = l-[(x + y)/At] [1 - exp(-At/x)] + (y/At)C 2c (x, y). 

In analogy with the expression proposed for the fluid velocity seen, cf. Eq. (|86|) . the following 
correction step is proposed, 
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For the simulation of the stochastic integral, one has {G 2 ,i is the A/"(0, 1) random variable used in 
the simulation of T" in the Euler scheme, see Table 0]) 
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where the second order moments ((f™ +1 ) 2 ) and (r™ +1 7™' +1 ) are computed from Eqs. (|141j) and 
()143p . respectively, by inserting the suitable timescales and diffusion matrix, that is r™, J'™+ 1 and 
B* . This completes the weak second order scheme. 

It can be shown, by means of stochastic Taylor expansion |17) . that the present scheme is a 
weak scheme of order 2 in time for system (|40|) . It is worth emphasising that no correction is done 
on position, x p (i), since the prediction is already of order 2. The complete scheme is summarised 
in Table El 

Limit cases: in limit case 1, when 1 <C At/r p , one has A2 C (t p , Ti) — > A 2 {At, Ti), B2 C (t p , Ti) — > 
£?2(Ai, Ti) and C2 C {j p , Ti) — » exp(— At/Ti). For the stochastic integral, one can show that — > 
7™ +1 . Inserting these results in Eq. (JS5J yields U^ 1 = U^f 1 , which is consistent with Eq. (|4i~j) . 
This result is a second order scheme for XJ p (t), and therefore the scheme remains of order 2 in 
limit case 1. 

When 1 <C At/Ti and Bi Ti = est (limit case 2), one has A 2c (r p ,T i ) — > A 2 (At,T p ) and 
B2c(T p ,Ti) — > B 2 (At,T p ), which gives for the numerical correction of the velocity of the parti- 
cles 

= \ U p,i exp(-At/r p ") + \ E£ 4 exp(-AVf; +1 ) 

+ A 2 (At, rZ)[(U?) + T % A?} + B 2 (At, f; +1 )[(U^ +1 ) + f p " +1 A« +1 ] ^ 

, -pn+l 
' i 

For the simulation of the stochastic integral, one can prove by looking at the limit values (when 
1 < At/T t and B % T t = est) in Eqs. (fT^Uj) . (IT4T|) and (TIlHl) that (here Q' pl is a W(0,1) random 
variable) 

f ^ J [ ' 2 ; n j [1 - exp(-2AVr»)] ^, (96) 

which is in accordance with Eq. ([72^1 . Unfortunately, it can be established, again by Taylor 
stochastic expansion, that the scheme is not of second order in time for system (|42l) , but of first 
order. This is due to the treatment of the correction step for the stochastic integral Fj(i) where 
r p has been retained in order to avoid anticipation and inconsistent numerical expressions of the 
Ito integral. As far as the fluid velocity seen is concerned, one has 




Kt 1 - <t/f +1 > + y ^" 2 ^ Qx,i, (97) 

which is in line with the previous result. This scheme is of second order, but the whole scheme is 
not. Indeed, as mentioned above, the scheme is only of first order for the velocity of the particles. 

When both the fluid velocity seen and the velocity of the particles become fast variables (limit 
case 3), that is when 1 <C At/Ti, 1 "C At/r p and Bi Ti = est, one can write for the velocity of the 
particle, for example from Eq. with 1 <§; At/r p , 



/ [R* T n+1 l 2 

u;T = (u- +1 ) + t; +1 a- +1 + x j y y w J g' Pti . (98) 



For the fluid velocity seen, Eq. (|97[) is unchanged. These results are consistent with the expressions 
of Section 15 .3. 31 In limit case 3, the numerical scheme for the position of the particles is equivalent 
to the Euler scheme written previously and is of first order in time. 

When the flow becomes laminar, that is when Ti — > with no condition on the product Bi Ti, 
one has the following limits: A2(At,Ti) —> 0, B2(At,Ti) —> 1 and %(t) — > 0, which gives for the 
fluid velocity seen, C/™^ 1 = (U" +1 ). For the velocity of the particles, the coefficients have the 
following limits: A 2c {t p ,T 1 ) -> A 2 (At,T p ), B 2c {T p ,Ti) B 2 (At,T p ) and C 2 c{T p ,Ti) which 
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gives together with the limit r,(i) — > 0, 



U n p f = 2 U v,i «p(-At/r^) + - C/ p " t exp(-Ai/r p " +1 ) 
+ A 2 (At, r;)[([/») + r" ^] + B 2 (At, f; +1 )[(U? +1 ) + r } 
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It can be shown, by regular Taylor expansion, that this scheme, together with the prediction step 
(Euler scheme) is a second order scheme for system (|4"4"|). 

In summary, a weak second-order scheme for system (1401) has been derived. This scheme 
satisfies all conditions listed in Section 15.31 However, second order convergence is not obtained 
in limit cases 2 and 3. In this latter case, the first order convergence is inherent to the spirit of 
the scheme, that is a single step to compute position x p (i) (in order to minimise the number of 
particle localisations in the algorithm). 

6 Specific and open issues 

The main objective of this paper is to present a consistent and rigorous numerical method for the 
computations of polydispersed turbulent two-phase flows using a mean- field/PDF approach. The 
mathematical framework and the models used in this approach have been put forward (Sections 
[2] and |3]) and, a general methodology, for the derivation of weak numerical schemes for the set 
of SDEs describing the dynamics of the stochastic particles, has been given in the context of 
particle-mesh methods (Sections H] and [5]). 

In this general methodology, the derivation of the weak numerical schemes has been performed 
only in the case of one-way coupling. As mentioned before, this is not a limitation of the method- 
ology and it is simply the status of the developments so far. The extension of the present results 
to two-way coupling is now discussed. Two issues are addressed: 

(i) the computations of the source terms in the PDEs describing the dynamics of the fluid mean 
fields, i.e. when two-way coupling is accounted for, that is when the particle mass fraction 
is high enough and the influence of the particles on the fluid mean fields must be taken into 
account, particle source terms are added to the fluid equations, cf. Eqs. (|129[) and ()130|> . 

(ii) the extension of the methodology introduced in Sections [4] and [5j i.e. when two-way coupling 
is considered, an acceleration is added to the SDE describing the dynamics of the fluid 
velocity seen, cf. Eqs. (|28]l and ([29)) . and the structure of the system of SDEs is changed. 

The first point is a specific issue, that is a practical solution is given for the computational proce- 
dure of the source terms. The second point is considered as an open issue since only explanations 
on the procedure to follow, for the extension of the weak numerical schemes to two-way coupling, 
are provided. 

After the treatment of the two-way coupling issues and before showing computational examples, 
some possible improvements and some open questions related to the numerical method will be 
discussed. The list of open questions related to the present numerical (particle-mesh) method is 
long. Here, attention is focused on two open issues: 

(iii) the formulation of boundary conditions in wall-bounded flows, 

(iv) the development of new numerical methods based on the present one. 
We start now with aspects related to two way coupling, i.e. issues (i) and (ii). 

6.1 Computation of the source terms 

As can be seen from the equations given in Table [TJ i.e. Eqs. (|129p and (|130l) . two source terms 
are present when two-way coupling is considered. The first one, say Sjj, represents the exchange 
of momentum between the discrete particles and the fluid. In the present paper, the only force 
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exerted by the fluid on the discrete particles is the drag force, see Eqs. ([T|), and by reaction the 
force exerted by the particles on the fluid is the reverse drag force. The mean momentum source 
term is expressed by 

(S u ) i = X ( Up>i ~ Us ' i ), (100) 

Tp 

where x = ( a pPp)/( a fPf)- From the particle equation of motion, cf. Eqs. (|29j) . the drag term 
is equal to the discrete particle acceleration (when gravity is first subtracted) and Su can be 
re-expressed as 

(Su)i = -X{^)- (101) 

From the discrete point of view, if we use the NGP technique, as explained in Section l5.2.1[ for the 
sake of simplicity (since most of what is presented below concern particle instantaneous quantities 
that can be put within the CIC formalism), this source term is the sum of the reverse drag force 
due to the discrete particles that are found in a given fluid cell. Then, the total fluid momentum 
in a cell [k], whose volume is can be written as 



N k (jjl \n+l _ (jjl \n 

fa \[k] ST^ l \ u p,i) \ u s,i) /ino\ 



where m p and U l p i stand for the mass and the velocity of the discrete particle labelled Z, respec- 
tively. The sum is performed over the N% particles that are located in cell [k] at iteration n 
(t = nAt). 

The source term for the fluid Reynolds-stress equations, see Eq. (|130p . raises new questions 
and its numerical evaluation is an interesting example of the specificities of stochastic calculus, 
cf. Section I2.2I For the discussion of its expression, we limit ourselves to the simplified case of a 
stationary one-dimensional system and to the source term, Sk, for the fluid kinetic energy k. The 
system of SDEs that we consider is 



dn p (t)= u '®- u >®*, 



u >tt) ^fU s (t)-U p (t) 



_ (103) 
dU s (t) = - — ^ dt-x [ ^ Sy " J T ~ PK "' j dt + VKdW{t), 

and the fluid kinetic source term which represents the work performed by the drag force is 

S k = x (U,(Z^-y (104) 

For this simplified case, and when the coefficients of the equations are constant, we can derive the 
analytical expression of the second-order moments. Indeed, for a long-enough time after the initial 
conditions, the stochastic process Z(t) = {U p (t),U s (t)} (since here x p (t) is irrelevant) reaches a 
stationary state and (Up), (U p U s ) and (C/f) become constant. Therefore, using Ito's calculus, we 
have that 

d{U 2 p ) =2(U p dU p ) =0, 

d (U p U s ) = (U p dU s + U s dU p ) = 0, (105) 
d (U 2 ) = 2 (U s dU s ) +Kdt = 0. 
The first two equations yield the equilibrium formulae for the second-order moments 



(U 2 ) = (U P U 
(UpU s ) = (U 



i 1 (1° 6 ) 
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while the third one gives the expression of the diffusion coefficient K to maintain a constant value 
of the fluid kinetic energy 

*- 2 <°?>(? + r£i)- (107) 

Using these formulae, the (equilibrium) analytical expression of the kinetic source term can be 
written as 

Sk = -x{U 2 s ) 7f ^—. (108) 

This source term is always negative which indicates that the drag force, which is indeed a friction 
force, induces a loss of energy in the fluid energy budget. This is valid for the total energy budget, 
whereas if we consider the fluid energy spectrum and its modulation by particles, particles may 
enhance turbulence at some lenghtscales (or wave numbers) due, for example, to wakes generated 
behind the particles. In the present model, we consider only Sk which is the integrated value of 
the exchange term over the whole spectrum, and, if we leave out the (possible) energy injected 
from particles by their initial conditions, the total energy gained by the particles comes from the 
fluid and the fluid kinetic energy source term is negative. However, when t p — > 0, that is when the 
discrete particles behave as fluid elements (but with a constant mass fraction, \), we expect the 
kinetic source term to vanish (Sk —> 0) since we consider a stationary case. Yet, from Eq. 1(1080 . 
it is seen that the limit is 

S k r>-x(U?)L (109) 

This spurious non-zero limit for vanishing particle characteristic timescale can be traced back to 
the Langevin model and is related to the fact that acceleration is indeed replaced by a white-noise 
term , cf. Chapter 6.8 in Ref. [8]. 

Nevertheless, it is possible to retrieve the correct limit in the numerical evaluation of Sk by 
resorting to a discretisation based on the Stratonovich definition, see Section 12.21 The first step 
is to write the source term with the particle acceleration as 

s k = -x(u.^). (no) 

Therefore, if we consider the integration of the source term in a time interval, we get 

S k dt=- X (U s dU p ) (111) 

and, in a formal sense, when t p — > 0, we expect the source term to become 

S k dt^ -x(U s dU s ), (112) 

since in that case U p — > U s , cf. Eq. (|4"Tj) in Section l4~2l Now, from Section [2~2l we know that 
the above expression can have different meanings. If we decide to regard the term (U s dU s ) as 
being defined in the Ito sense, as it should be in order to be consistent with the algebra retained 
throughout the paper, we would find the non-zero limit given above. Yet, if for this expression of 
the source term, we decide to consider it as being defined in the Stratonovich sense, then 

S k dt ► -x (Us o dU s ) = - J (d(U s f), (113) 

which is indeed zero since we are in a stationary case. 

The difference between the two stochastic calculus is only presented here since it provides a 
useful guideline. In the present case, it is seen that the interest of the Stratonovich expression is 
that the formal quantity dU s /dt can still be handled as if it were a normal derivative (and, in our 
case, the limit of dU p /dt when r p — > 0). This suggests therefore to express the kinetic source term 
numerically, in a fluid cell [k] at time t = n At with Nk particles, as 

N k 1 _ (Tjl\n 

a fPi vf S« = - £ ml 1 - ((Ulr +1 + Mr) ^ ^ ■ (114) 



37 



From the properties of the numerical schemes developed in the previous Sections, we have that 
{U l p ) n -> {U\) n when t p -> 0. Thus, in that limit 



N >= - 'rjl\ n +l _ (Tjl\n 



^PfV l ; ] si k] ^-J2m l p U(U l s) n+1 + (Ui) 
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1=1 

and, when the stationarity of U s is indeed enforced numerically, this term is zero. Finally, going 
back to the exact fluid Reynolds stress equations, we propose to express the numerical source 
terms as 

N k ( . rjjl \n+l _ (Tjl \n 
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As mentioned at the beginning of this section, most of the arguments that have been presented 
concern the discrete evaluation of each particle term U s (dU p /dt). The expression proposed above 
in the NGP formulation can be used directly within the CIC technique. Another interesting 
question is to ask to what cell (or cells) the different source terms should be assigned. Indeed, 
within one time step, particles may cross several fluid cells and the source terms, say Sjj and Sr 
which represent the total momentum and energy exchange terms, should be distributed between 
the different cells crossed by the particles. One possibility is to apply to each fluid cell crossed 
by a particle, the different reverse expressions, say dU p /dt and U s (dU p /dt) in proportion of the 
time spent in that cell (the residence time) which is then a fraction of the time step. This is 
probably the most precise expression and the most accurate discrete formulation, but it implies 
to keep track of the different fluid cells along the particle trajectory within one time step. In 
a complex geometry and unstructured meshes, given present localisation algorithms, this is not 
an easy task and it induces computational overloads. For these reasons, at the moment, it is 
proposed to evaluate the total source terms from the particles that were located in that cell at 
the beginning of the time step. This evaluation has been applied in the various computational 
examples presented in Section [7] It can be seen as a first-order spatial approximation or based on 
an implicit assumption that the particle Courant number remains of order one in most cases. 



6.2 Extension of the weak numerical schemes 

When two-way coupling is accounted for, the SDE describing the dynamics of the fluid velocity 
seen is supplemented with an acceleration term, cf. Eqs. (|28|) and ([29)) . in order to account 
for the influence of the discrete particles on the statistics of the fluid velocity sampled along the 
trajectory of a discrete particle. This supplementary acceleration changes drastically the nature 
of the equation system and one has (the equation for position is omitted for the sake of clarity) 

( 11 

dUp.i(t) = U p .i dt H U s ,i dt + gi dt, 

T p T p 

X / l v \ ^ (117) 

dU s>i (t) = -Up,, eft - — + — U Sli dt + Cidt + Y, Bij dWj{t), 

Tp \li Tp J 

V J 

that is, the SDE for the fluid velocity seen, U s (t), depends explicitly on the velocity of the discrete 
particle, U p (t). This dependence complicates the analysis of the system, in particular the limit 
systems when the time scales go to zero, cf. Section [4] If one is able to find the limit systems 
in the continuous sense, the extension of the numerical schemes can be obtained in the same way 
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as presented in Section 15.31 However, in order to calculate the analytical solution with constant 
coefficients, one has to express the following matrix in diagonal or triangular form (depending on 
the roots of the characteristic polynomial) 



-1/r, 
x/r P 



p 



(l/T t + X /r p ) 



(118) 



Once this is done, the previous analysis can be used, but in the frame of much more complex alge- 
bra. Once the analytical solution is obtained for the eigensystem, one has to go back to the original 
system (state vector) with some transformation matrix (which is formed by the eigenvectors). 

This difficulty is actually not a typical feature of two-way coupling. As a matter of fact, in 
the one-way coupling case, if the alternative model is chosen, cf. Section 13.1. 2\ the drift term is 
written in terms of the local instantaneous velocities. Therefore, an acceleration which has the 
same form introduced 



where the time scales of the mean flow, Ty, are given by T-" 1 = d(Ui}/dxj. As a consequence, 
if such a model is used for the drift term, the problems inherent to the derivation of the present 
schemes with two-way coupling are already encountered for the one-way coupling case. There is, 
to our knowledge, no specific work in the literature dealing with this subject. 

6.3 Boundary conditions in wall-bounded flows 

In the present work, for the computational examples, cf. Section [71 the wall boundary conditions 
for the system of SDEs are treated as follows: for the discrete particle velocity, TJ p (t), an elastic 
wall-particle collision is applied whereas for the fluid velocity seen, U s (t) , we build on ideas from 
turbulent single-phase flows i36j 49 in order to ensure consistency when t p — > 0. As a matter of 
fact, in one-point PDF models for single-phase turbulent flows, cf. Section 13.1.11 or in one-point 
PDF models for the discrete particles, cf. Section l3.1.2[ the derivation of boundary conditions for 
fluid or discrete particles, when solid boundaries are present, has not received the needed attention. 
Here, for the sake of simplicity, we make our point by considering, as an example, only the motion 
of fluid particles. 

In the framework of PDF methods for single-phase turbulent flows, boundary conditions of the 
wall function kind have been studied and proposed 36, 50 . This solution has been investigated 
rigorously from the mathematical and physical (based on the knowledge of the phenomenology of 
the near-wall region) points of view. In practice, the numerical treatment is developed in analogy 
with the wall-function approach used in RANS computations, a method which is perfectly in line 
with one-point high-Reynolds number PDF models. However, in some engineering applications 
where a precise description of the near-wall region is needed, it may be of interest to replace the 
wall- function boundary conditions with a direct particle- wall interaction, i.e. U/(i) = at the 
wall. According to Section 13.1. 1[ a one-point PDF model for single-phase turbulent flows reads 



and it has been shown in Section 15.3.41 that a possible weak second-order scheme is (when *Bjj 
verifies property ([75)) ) 
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where = %,{t + Ai, J/™/ 1 ) and Q3^ +1 = T> i:j (t + Ai,^ 1 ). This scheme is used for 

our present discussion and it is different from the one developped above, cf. Table [SJ With this 
stochastic framework in mind, some open questions remain. 

(i) Is it possible to propose a general form of wall boundary conditions for the fluid particles, 
independently of our particular Langevin model? 

(ii) What boundary condition ensures that the impermeability condition is valid just at the wall 
(as in the real world)? 

(iii) What is the order of accuracy of the boundary conditions in the frame of our numerical 
schemes? 

(iv) Is it possible to propose high-order (second-order) boundary conditions? 

The above questions might seem easy to answer at first glance, but the subtleties of stochastic 
calculus make these open issues difficult to solve. At present, only one proposition has been 
made [51] . In that work, the authors have proposed to impose a zero velocity to the particles 
reaching the wall during a time step and to move them in space by symmetry at the wall. The 
proposed treatment is sensible, but it presents some shortcomings, i.e. it remains dependent on 
the particular model used (a Wiener process was used in the equation for position to account for 
viscous effects in the vicinity of the wall), and the order of accuracy is not given. A mathematical 
approach to this problem can be found in the book of Ottinger |31j . 

The scarcity of the literature on this subject calls for future rigorous development in the 
formulation of boundary conditions at the wall in one-point PDF methods. A good illustration 
of the lack of knowledge will be presented in Section [7] for a computational example of particle 
deposition phenomena. 

6.4 New hybrid methods 

In the present work, a hybrid method has been used: the fluid is described with a mean-field 
(RANS) approach whereas the statistics (the pdf) of the discrete particles are reproduced by 
introducing stochastic particles (SDEs). 

In stand-alone methods for one-point PDF models for single-phase turbulent flows, it is known 
that the bias (cf . Section I5.1.2j) is the main concern in the control of the numerical error [41] . If 
this is also the case for the numerics put forward for one-point PDF models for discrete particles, 
every idea improving this shortcoming is welcome. A solution could be to resort to VRT that 
have been developed in disparate fields. A possibility could be to resort to a hybrid algorithm for 
the numerical treatment of the discrete particles where some variables could be solved by a mean- 
field method and others by a PDF method. In such configurations, duplicate fields usually arise, 
and consistency conditions must be imposed. These consistency conditions give the opportunity 
of introducing VRT. Indeed, the mean variables computed from a mean-field method are by 
construction not biased. If the PDF method contains the evolution in time of the corresponding 
instantaneous variables, the operation of centering the moments extracted from the PDF approach 
with the ones computed from the mean-field algorithm leads to an excellent reduction of variance 
|37j . Moreover, it would be helpful to find a criterion, in the frame of domain decomposition, 
in order to use the mean-field or PDF algorithms where it is most appropriate. For example, in 
some parts of the flow where the knowledge of some mean-fields is sufficient, one would resort 
to mean-field algorithms whereas in other regions, where the physics are complex and the pdf 
is needed, one would use a PDF algorithm. In such an approach, the central issue becomes the 
consistency at the boundaries between the contiguous domains. Some work in that sense has been 
carried out in the field of Direct Monte Carlo Simulations (DSMC) 52 . 

At last, in the present work, the proposed particle algorithm (the numerical algorithm for the 
set of SDEs, cf. Section Rj.l.ip is compatible with other approaches for the fluid. It is one of the 
strong points of this numerical method. Therefore, it is conceptually possible to think about some 
other configuration and in particular to a LES/PDF one, i.e. the set of SDEs is provided with 
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filtered fluid fields instead of mean fluid fields. Even though, in such a configuration, an increase 
of the computational effort is expected, the quality of the results, in particular, for cases where 
RANS models are known to be inadequate, could be improved, cf. the computational example 
for particle deposition in Section [7J In such an algorithm, the challenge is to reconstruct the 
subgrid scale fluid velocity along the discrete particle trajectories [6j [T] [53] . A possibility is to use 
PDF methods; it has been attempted in single-phase flows [54] and it remains to be developed for 
dispersed two-phase flows. 

7 Computational examples 

Three numerical computations of polydispersed turbulent two-phase flows are now presented. The 
first one (swirling flow) is chosen in order to show that significant improvements in the computing 
efficiency can be achieved by using, for the integration of the set of SDEs, a second-order scheme 
instead of a first-order scheme. The second computation (bluff-body flow) demonstrates the ability 
of the models to capture the main physics of the flow and the specificity of PDF models from 
which valuable information can be extracted. The third example (particle deposition) illustrates 
the ability of PDF models to treat flows were complex physics are involved. 

In the first and third examples, i.e. swirling flow and particle deposition, both flows are dilute 
enough so that only one-way coupling is under consideration. The numerical schemes presented in 
Section [5. 31 can therefore be used directly. In the second example, bluff-body flow, the suspension 
is rather dense and one has to take into account two way-coupling. Numerically, for the integration 
in time of the set of SDEs, this is done by resorting to the first-order scheme and by treating the 
coupling term, cf. Eq. (|28|) . as an explicit source term. Collisions, that might occur in some 
restricted areas of the computational domain, are not taken care of. It is, however, fully possible 
to treat collisions between discrete particles in the frame of the present PDF approach. This 
has been discussed elsewhere [8] on a theoretical basis and the inherent numerical developments 
remain to be done. 

7.1 Swirling flow 

In this particular example, no comparison with experimental data is attempted since the purpose 
of the computations is to show the benefits of using second-order schemes instead of first-order 
schemes for the integration of the set of SDEs. 

7.1.1 Experimental setup 

The turbulent polydispersed two-phase flow under investigation corresponds to a gas-solid flow 
(air and solid particles) in a cyclone of the Stairmand type [55], see Fig. [TO] Cyclone separators 
are devices used to separate particles from gas flows. The gas flow inside the cyclone has quite 
complicated patterns, that is a reverse swirling flow with quite high rotational velocities. The swirl 
is created by the tangential inlet, but it is well-known from experiments, that the gas flow exhibits 
a double helix structure. The flow spirals downwards (with a constant intensity) to the vortex 
finder (exit tube of the cyclone at the bottom) where it reverses and spirals upwards in a cylindrical 
volume having roughly the diameter of the exit. In such a device, the separation between air and 
particles is not due to gravity but to the effect of the double helix. Indeed particles entering the 
device are entrained towards the outer wall (by centrifugal forces) where they flow downwards to 
the exit (the axial velocity of the gas is oriented downwards at the walls) . 

The efficiency of a cyclone is characterised by its selectivity curve. This curve expresses the 
ratio (in mass) of captured particles as a function of their diameter. For very small particles, this 
curve goes to zero efficiency as their inertia decreases (t p — > 0), i.e. particles tend to behave as 
fluid elements. On the contrary, large particles are all collected and the efficiency converges to 1. 
Between these two asymptotic cases, the cyclone efficiency is an increasing function of the particle 
diameter. 
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7.1.2 Numerical simulations 

The simulated cyclone [SHI [S7] has a diameter D = 0.2 m, the glass particles have a density 
of p p = 2500 kg/m 3 and diameters ranging between 0.5 and 5 /im. The gas (air at ambient 
temperature ~ 293 K) is injected with a constant velocity of 30 m • s -1 . In the present simulations, 
the flow is dilute enough (the mass of particles per unit of gas is quite low) not to consider 
two-way coupling effects (the particles have no influence on the flow field). In addition, it is 
well established from experiments that such flows are stationary. Consequently, the flow field is 
computed in advance and all particles are tracked in a frozen field. The prediction of the flow 
field is rather challenging given the complex structure of the flow. In this work, a second-order 
turbulence model (Rotta model [13], which is consistent with the form of the SDEs) was used 
with a fine grid (approximately 4 • 10 5 nodes) in order to obtain mesh-independent calculations. 
Figures [IT] and [12] show the axial and radial mean velocity profiles of the air flow at two different 
heights: it can be seen that the numerical results are in good agreement with the measurements. 

Particles are then tracked in this frozen field. The diameter range of the glass particles has 
been discretised as follows, d p — [0.5; 1; 1.5; 2; 3; 4; 5] /im. For each class (diameter of particles), 
a number N pc of particles is released (this number is identical for each class). The computation 
stops when all particles have left the computational domain. 

7.1.3 Results and discussion 

A numerical study has been carried out to show that the results (the obtained selectivity curves) 
are independent of the time step, At, and the number of particles per class, N pc . For instance, 
for the second-order scheme, the computations show, for two different time steps, that roughly 
400 particles per class are necessary to obtain selectivity curves which do not depend on N pc , see 
Fig. [T3J The time step of At — 10~ 4 s guarantees that the results do not depend on the time 
discretisation, Fig. [JjJ] For the first-order scheme, similar results are obtained : At = 5 ■ 10~ 6 s 
and N pc = 400 ensure that the results depend neither on the time step nor on the number of 
particles per class, Fig. [131 

It is then observed that, for this particular flow, there is a great difference between the re- 
spective time steps for the first and second order schemes, see Fig. I13[ and this for the same 
numerical results (selectivity curve). In fact, with a second-order scheme, the time step can be 
multiplied by a factor 20 compared to a first-order scheme. If one accounts for the computer time 
(the computation of a time step takes approximately 30 % extra time compared to the first-order 
scheme), there is a gain in CPU time by a factor 15 by using a second-order scheme. Therefore, in 
this case, it is seen that the complexity of the second-order scheme is balanced by the reduction 
of computing time. 

It can be stated that, for flows where velocities and the curvature of the trajectories of the 
discrete particles are important, it is recommended to use a second-order scheme rather than a 
first-order one, unless one is ready to pay the computational price. This is clearly seen in this 
computational example where the time step of the first-order scheme is extremely small: as a 
matter of fact, in such a flow, a precise prediction of the particle velocity is needed because the 
numerical error on this quantity amounts at simulating an additional centrifugal force. 

7.2 Bluff body flow 

In this second example, the numerical results are compared to experimental data in order to 
show the ability of the present approach to reproduce the main trends of complex turbulent 
polydispersed two phase-flows. Other results are displayed to enhance the specificity of the mean- 
field/PDF approach, that is the type of information which can be extracted. 

7.2.1 Experimental setup 

The 'Hercule' experimental setup [581 159] is characteristic of pulverised coal combustion where 
primary air and coal are injected in the centre and secondary air is introduced on the periphery, Fig. 
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[T4l This is a typical bluff-body flow where the gas (air at ambient temperature and atmospheric 
pressure) is injected both in the inner region (jet) and the outer cylinders (exterior). The ratio 
between the gas velocity in the inner region, Uj, and the gas velocity in the outer region, U e , is 
low enough so that a recirculation zone downstream of the injection is created. Two honeycombs 
are used in the experimental setup in order to stabilise the flow so that no swirl is present. Solid 
particles (glass spheres) are then injected from the inner cylinder with a given mass flow rate. 
The injected glass spheres have a density p p = 2470 kg/m 3 and a known diameter distribution, 
typically between d p = 20 /im and d p — 110 /im around a mass- weighted average of d p ~ 65 /im. A 
polydispersed turbulent two-phase flow which is stationary and axi- symmetric is then obtained. 
Moreover, two-way coupling takes place since the particle mass loading, (j> = x(Up)/(Uf), at 
the inlet is high enough. Experimental data is available for radial profiles of different statistical 
quantities at five axial distances downstream of the injection, Fig. [H] (axial profiles along the axis 
of symmetry have also been measured). The 'Hercule' experimental setup is a very interesting test 
case for polydispersed turbulent two-phase flow modelling and numerical simulations since most 
of the different aspects encountered in such flows are present. The particles are dispersed by the 
turbulent flow but in return modify it. Furthermore, the existence of a recirculation zone (with two 
stagnation points, Si and S2 in Fig. I14|) where particles interact with negative axial fluid velocities 
constitutes a much more stringent test case compared to cases where the fluid and the particle 
mean velocities arc of the same sign (the problem is then mostly confined to radial dispersion 
issues). These features are displayed in Fig. [Til where mean streamlines are shown (solid lines for 
the fluid and dashed lines for the particles). For the fluid, there is a rather large recirculation zone 
with stagnation points. For the particles, depending on their inertia, several behaviours can be 
observed: some particles do not 'feel' the recirculation zone and leave the test section immediately. 
Others are partially influenced and change direction before leaving the apparatus, whereas some 
particles follow closely the recirculation pattern. This will be seen in the results showing the pdf 
of the particle residence time at different locations in the flow. 

7.2.2 Numerical simulation 

A two-dimensional, single block, non-Cartesian, non- uniform mesh (142 x 3 x 75 nodes for the 
(x, r, 0) coordinate system) has been generated in accordance with the axi-symmetric property 
of the flow. It was carefully checked that the results are not too sensitive either to the time 
and spatial discretisations or to the number of stochastic particles (At = 10~ 3 s and N = 14000 
particles). A second-order turbulence model (Rotta model, which is consistent with the form of 
the SDEs) was used. The projection and averaging operators were approximated with a NGP 
technique [33]. For further details on the numerical computations, see Ref. [60 . 

In the simulations, the following procedure is adopted. The single-phase flow case is first 
computed until the stationary state is reached. By doing so, it is possible to check that the 
prediction of the flow field, without the particles, is accurate. Then the discrete particles are 
introduced until the stationary state is obtained again. At that point, the number of particles in 
the flow is roughly constant (it fluctuates around a mean value). From this state, computations 
are continued to extract the statistics which are compared to the experimental values. The last 
computation is performed to allow time averaging on the ensemble averages so that the statistical 
noise can be reduced to a minimum (VRT). 

7.2.3 Results and discussion 

Figures [T5] fT6l and [T7l show that, in this particular flow, there is almost no difference between the 
predictions with the first and second-order schemes. The main differences take place in regions 
where only few particles are present (large diameters) and consequently the statistical results 
contain some noise (see the ragged behaviour of the curves). 

Three sets of results are given: (i) radial profiles of the particle axial velocity, Fig. [15] (ii) 
radial profiles of the particle radial velocity, Fig. [16] and (iii) radial profiles of the fluid axial 
velocity, Fig. 1171 All sets of numerical results compare relatively well with the experimental data, 
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in terms of the shape of the curve and of the magnitude which are observed. The recirculation 
zone (x — 0.16 m) is well predicted as indicated by the velocity profiles. All (mean and fluctuating) 
velocities go to zero when no particles are present. The widths of the numerical curves indicate 
that the predicted radial dispersion of the particles is also in line with experimental findings. 

In Figs. [15] to [TTl only first and second-order moments of the variables of interest have been 
displayed (mean and fluctuating velocities for the fluid and the particles) . This information could 
have been obtained by resorting to classical mean-field equations. However, in many engineering 
applications (for example combustion), it is necessary to know the distribution of the residence 
time of the particles at a certain time. In other words, one would like to know for all the particles 
found in a certain zone how much time they have spent inside the domain, or even if they have 
previously entered a marked region. This kind of information is not available in a mean-field model 
whereas in the present hybrid approach, this information is directly provided without additional 
costs: the pdf of the variables attached to each particle, which contains far more information than 
a few moments, is explicitly computed. 

A typical example of this type of information is given in Fig. 1181 In the first plot, on the 
left-hand side, a snap shot of the local instantaneous positions of the particles is given, where 
particles are coloured by their residence time. Two distributions are extracted, one in a cell close 
to the inlet and the other one in a cell close to the outlet. The pdf in the cell near the inlet clearly 
shows the recirculation pattern: the distribution is highly peaked, which represents particles which 
have just entered the domain, but a small number of particles have a quite high residence time, 
i.e. they recirculate. At the outlet, since the particles have different trajectories in the domain, 
a continuous spread in residence time is observed. In the region near the inlet, more information 
can be gathered, for example the local instantaneous axial velocity, see the RHS graph in Fig. 
IT51 Most particles have the same axial velocity, around 4m/s which is actually the inlet velocity. 
These particles correspond to the peak observed in the residence time: they have just entered the 
domain and travel directly to this region. The rest of the particles have a smaller axial velocity 
but with much wider fluctuations: these particles correspond to particles which are recirculating. 

At last, it is often argued that the mean-field/PDF approach is time-consuming. As a matter 
of fact, in this computation, the time spent by both solvers (the mean- field solver for the fluid 
and the PDF solver for the particles) has been compared for the same number of computational 
elements (mesh points and particles, respectively). It is found that the PDF solver is slightly faster: 
this is not really surprising since the mean-field fluid solver implies the use of a full second-order 
turbulence model (6 coupled PDEs). 

7.3 Pipe flow: deposition 

In this last example of numerical applications with the mean-field/PDF approach, a flow where 
complex physics are involved, i.e. particle deposition, is under investigation. Particle deposition 
from a turbulent flow on walls is a phenomenon which is observed in many engineering applications 
(for example thermal and nuclear systems, cyclone separators, spray cooling) and also in various 
environmental situations. Given the large number of possible applications, a lot of interest has 
been devoted to this subject and many studies have been carried out in the last decades. 

7.3.1 Experimental setup 

Different experiments have been conducted to observe deposition in turbulent flows. In most of 
them, attention is focused on the deposition velocity |611 162) which is defined as k p — m v j C, 
where m p is the mass flux and C is the bulk mean particle concentration. This deposition rate, 
often presented as the dimensionless deposition velocity k p /u* , is a function of the dimensionless 
particle relaxation time, r+ defined as 
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where S + is the dimensionless stopping distance, U p q is the particle initial velocity and u* is the 
friction velocity, u* is evaluated with the Blasius formula, u* — [0.03955 Re ' 25 ] 5 U m , where U m is 
the bulk mean velocity. The deposition velocity is the key point in many engineering applications 
where one seeks the law that gives k p /u* as a function of r p , that is as a function of the particle 
diameter. Recently, several experimental studies and DNS studies of particle deposition have been 
presented, for example |63[ 164] , and have improved the understanding of the physical mechanisms 
at play. In particular, a lot of information has been obtained on the dynamical structures of wall- 
bounded flows, like the coherent structures which manifest themselves in the near-wall region. 

In the present computational example, the principal interest is to show the advantage of solving 
the set of SDEs with a numerical scheme consistent with all limit cases, see Sections l4l 15.3.31 and 
15.3.41 Indeed, in pipe flows with particle diameters ranging from 1 /im to 100 /im, all limit cases 
can be encountered: 

1. for the smallest particles (r p — > in the continuous sense and r p <C At in the discrete sense), 
limit case 1 is obtained, 

2. in the near wall region, i.e. T£ i <C At, for example in the peak-production region where 
turbulent kinetic energy is maximal, one has By T£ , = est, a situation which is characteristic 
of limit case 2, 

3. in the same region as above with small particles, one has limit case 3, 

4. in the close vicinity of the wall, i.e. T£ i <C At, with no condition on the other moments, 
the flow is laminar, that is limit case 4. 

The need to cope with limit cases is the result of practical considerations. Indeed, if the nu- 
merical scheme were not consistent with the limit cases, it would be very inefficient to carry out 
computations with a time step limited by the smallest timescale. 

In the present work, numerical simulations corresponding to the experimental setup of Liu and 
Agarwal [5T] are presented, i.e. the deposition of particles (920kg/m 3 in density and diameters in 
the range 1.4 to 68.5 /um) in a vertical pipe flow at a Reynolds number of 10 4 . 

7.3.2 Numerical simulation 

In order to describe the particle phase, 10 4 stochastic particles (distributed in 10 diameter classes, 
cf. Table O are released in a frozen field, i.e. the flow is stationary and dilute enough so that only 
one-way coupling is under consideration. Two frozen fields are computed with a standard k — e 
turbulence model and with a Rij — e (Rotta model) both with wall-function boundary conditions. 
The computations are performed with a 2D mesh, 100 x 20 x 3 nodes (the flow is axisymmetric) . 

To compute the deposition velocity, the fraction of particles remaining in the flow, F, is 
evaluated as a function of the axial position x [65j . F is calculated as the number of particles that 
reach the sampling cross-section, divided by the total number of released particles. The particle 
deposition velocity is then computed as follows [55] 



where Dh is the diameter of the pipe and Fi is the value of F at a given sampling cross section 
labelled i (axial position xi). 

Numerical tests have been performed to check that the numerical results are independent of the 
values of the numerical parameters, in particular the number of particles, N , and the time-step, 
At. It was checked beforehand that the numerical prediction of the fluid field is grid-independent. 

Both numerical schemes (first and second order) were tested with different time steps, cf. 
Fig. 1191 All computations were then performed with the weak second-order scheme and a time 
step of 10~ 4 s. Indeed, Fig. [T9l shows that both schemes give similar results and that a time 
step At — 10~ 4 s ensures that the computations are independent of At. The independence of 
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the deposition velocity for the whole range of particle diameters with respect to the time step 
illustrates the benefits of a numerical scheme which is consistent in all limit cases: for instance, 
the values of the particle relaxation timescales given in Table [5] cover three orders of magnitude 
(limit case 1). Computations can anyway be carried out with by using the same constant time-step 
for all classes and in the whole domain. 

An analysis of the statistical error has also been carried out. Since particle deposition velocities 
are calculated by a Monte Carlo method, it is important to check that the number of particles 
(which represents samples of the pdf) is sufficiently large so that statistical error is reasonably 
small. In Fig. [201 results obtained with three different values of N pc (the number of particles used 
for each class of diameter) are presented (N pc = 500, 1000 and 5000). Fig. [20l shows that there is 
no clear difference between the results obtained with different values of N pc . As a matter of fact, 
it seems that 500 particles for each class of diameter is enough. Nevertheless, for all following 
simulations, the value of Npc = 1000 particles for each class of diameter has been chosen, in order 
to reduce statistical noise. 

7.3.3 Results and discussion 

The numerical results are now compared to experiments and some sensitivity tests are conducted. 
Some proposals are put forward for the features which seem to be the most significant ones for 
a good representation of deposition phenomena. Possible improvements of the computational 
method shall be exposed. 

In Fig. [21] results obtained with the general PDF model, Eqs. (dJ and ([24]) . and two turbulence 
models (k — e and Rij — e), are displayed. The difference between the simulations performed with 
the two different turbulence models is negligible: this is not too surprising, since for turbulent 
pipe flow, both models give similar mean fluid velocity profiles. The standard PDF model is 
integrated with wall-function conditions for the fluid and pure-deposition boundary conditions for 
the particles. These results are coherent with those obtained in an analogous configuration by 
Schuen 166 . 




Figure I2T1 shows that, for heavy particles (r+ > 10), the model predictions are in good agree- 
ment with experiments whereas for light particles (r+ < 10), the deposition velocities are strongly 
overestimated (they remain at the same level as that of the heavy particles). Therefore, the model 
is not suitable for simulations of deposition phenomena in the range r+ < 10. This statement 
is consistent with experimental and DNS findings [67]: heavy particles are slightly affected by 
near-wall boundary layer and more especially by the specific features of the local instantaneous 
turbulent structures in the near-wall region. On the contrary, for light particles, the physical 
mechanism of deposition changes, with a growing importance of turbulent structures and near- 
wall physics. In the current PDF model, near wall physics are mainly described by wall-function 
boundary conditions which may be sufficient for heavy particles deposition but not for light particle 
deposition. 

Wall-functions give a reasonable approximation of the mean fluid velocity profile in the loga- 
rithmic region, but in any case they do not describe the viscous sub-layer. Therefore, a question 
arises: is the prediction of small particle deposition velocity sensitive to changes in the fluid mean- 
field profiles? This matter was recently investigated for other Lagrangian models |65j . Following 
the same reasoning, simulations have been carried out with a given frozen field (consequently, 
wall- function boundary conditions are suppressed). The frozen field can be obtained, in this par- 
ticular case, either from analytical solutions for the mean fluid fields ((U), k, (e)) [4] and/or from 
DNS data. In Fig. [22] two frozen fields are tested. In the first field, the axial mean fluid velocity 
(Ui) is given by the law-of-the-wall equations (the values k and e are that of the computations). 
In the second field, (Ui) is still given by the law-of-the-wall, and the turbulent kinetic energy, k, 
and the turbulent dissipation rate, e, are curve-fitted to the DNS data that can be found in the 
work of Matida et al. 65 . Thus, in the second field, mean fluid profiles are exact. Figure [22] 
shows that an exact frozen field hardly improves the results. An explanation might be that the 
eventual effect of the exact mean fluid profiles are concentrated in a very thin region. The most 
important quantity is expected to be the turbulent kinetic energy, which goes to zero at the wall 
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and should affect mainly light particles. Nevertheless, k diminishes only from y + sa 10, where it 
has its maximum (peak production). The resulting effect is not easy to be foreseen and it may be 
negligible with respect to the overall effect of migration of particles towards the wall due to the 
net mean flow. 

In order to further support the argument above, mean near-wall residence times of deposited 
particles have been computed in the layer y + < 30, for each class of diameters. Indeed, this 
quantity has been found to properly distinguish different deposition mechanisms |67j . A rough 
description of the physics of deposition is that heavy particles, which are slightly influenced by near- 
wall structures, deposit with small near-wall residence times by the so called free-flight mechanism. 
On the contrary, light particles are trapped and driven by turbulent structures and depose with 
large near-wall residence times, this mechanism is called diffusional. The lighter the particles are, 
the more important the diffusional mechanism is. In Table |5J the results obtained for each class of 
diameters are given for a simulation corresponding to the exact frozen field. For the sake of clarity, 
the residence time is always expressed in non-dimensional form (it is normalised with the viscous 
timescale, Vf/(u*) 2 . t + = t{u*) 2 /vf). Table [6J shows that all particles deposit after small near- 
wall residence time, that is by free-flight mechanism. Moreover, since the residence time slightly 
increases with d p , the motion of particles is influenced by the migratory flux. The force exerted on 
the particles by the fluid being inversely proportional to d p , light particles reach the walls faster 
than the coarse ones. Therefore, in absence of a representation of turbulent coherent structures 
(which should be able to trap particles in the near-wall region and which should describe correctly 
the mechanisms of deposition) , the sole mean fluid profiles are not the main mechanism. 

Two possibilities exist to improve the prediction of deposition phenomena and in particular of 
light particles: 

(i) some phenomenological model can be introduced, based on the present knowledge of de- 
position physics. In this type of approach some hypotheses are made in accordance with 
experimental findings. Some parameters may be present and may be fitted in order to find 
good comparison with experiments. This approach ought to verify if the hypotheses made 
are correct or not, and, thus, ought to show which are the dominant aspects not covered by 
the standard model. 

(ii) it is also possible to propose extensions of the present PDF model that reproduce correctly 
the variations of the fluid statistical moments (such as (U), Rij, e, etc.) throughout the near- 
wall region, including the viscous sublayer |49) . Such a model might lead to improvement 
for the deposition velocity. However, it would require a very refined mesh in the near-wall 
region, considering the high value of the Reynolds number. Furthermore, only the statistics 
of the fluid would be well reproduced and, as mentioned above, it is believed that the 
contribution of the specific features of coherent structures should be considered for small 
particle deposition. 

Therefore, for practical purposes, the first proposition has been retained [55] ■ This subject is not 
further developed here and it is left as a challenging open issue. 

8 Conclusions and perspectives 

In this paper, we have presented a comprehensive review of the numerical methods involved for 
the computation of polydispersed turbulent two-phase flows using a particle stochastic method 
based on Langevin equations. The present mean-field/PDF model is one among a host of other 
so-called Euler/Lagrange models. However, it is worth putting forward two main specific aspects 
of the present framework. 

(i) The usual term Euler/Lagrange refers to the point of view adopted for the description of 
the two phases: an Eulerian point of view for the fluid phase and a Lagrangian one for the 
particle phase. This terminology can be misleading and does not clearly identify the physics 
involved, the level of information contained in the statistical description and the numerical 
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tools which are adopted. For example, the so-called Eulerian equations for the fluid phase 
can be directly obtained through the two-particle stochastic formulation sketched in Section 
13.2. f l and Appendix^ and its numerical solution may involve Lagrangian ideas (for instance, 
in the method of characteristics for the discretisation of the convection terms). And, at least 
in theory, the PDF equation could also be solved using a mesh and an Eulerian description in 
phase space. In the present work, the complete model is called a mean-field/PDF approach. 
This refers directly to the level of information contained in the description: the fluid phase 
is described by a limited number of statistical moments (in practice, at most two for the 
fluid velocity) while the particle phase is characterised by the PDF of the variables retained 
for its description. The complete model is therefore a hybrid model. The hybrid nature of 
the model is then reflected in the numerical approach developed in Section The fluid 
mean fields are computed as the solution of PDEs, involving a mesh and, say, classical 
Finite Volume schemes. As explained in Sectional the PDF equation is solved, in a weak 
sense, by a particle Monte Carlo method where the particles should be seen as instantaneous 
realizations of this PDF rather than real particles. The numerical approach is thus an hybrid 
PDE/Particle Stochastic method, or a Mesh/Particle Stochastic method; its specificities have 
been discussed at length in Section [5] The various details treated in that Section can be 
improved, but the important point is that they are developed once a clear framework about 
the complete hybrid method is first set forth. This framework is a necessary guideline. 

(ii) Drawing on these first remarks, the second aspect of the present formulation is the fact that 
there is a separation between the theoretical construction of the model and its numerical 
solution. Indeed, given the correspondence (in a weak sense) between the PDF equation 
and particle stochastic equations, as explained in Section [51 it may be tempting to develop 
directly the model in discrete time. This amounts to treating without distinction the model 
and its numerical scheme. This may be confusing and may not help to identify the actual 
issues. In the present work, the theoretical stochastic model is developed and first written 
in continuous time. This requires knowledge of the mathematical background of stochastic 
diffusion, but actually this effort is a valuable investment and the formulation in terms of 
the particle trajectories of the stochastic process in continuous time simplifies the situation 
and the numerical developments. 

The central point of the complete work presented in Section[3]concerns the multiscale character 
of the stochastic theoretical model which is then reflected in the numerical developments of Section 
O Placed between the mathematical background and the numerical implementations, this Section 
illustrates the interplay between mathematical formulation, physical modelling and numerical 
developments. The mathematical manipulation of the system of equations reveals the property of 
the model: different limits are continuously reached depending upon the values of the observation 
timescale with respect to the different physical timescales. These limits correspond to natural 
physical diffusive limits and they point out that some physical variables are not real white-noise 
terms, but that their effects may be regarded as such at a certain scale. Using the numerical 
time step as the observation timescale, this physical property appears in turn as a basis for the 
development of the numerical schemes which, while being explicit and stable, can satisfy these 
limits without any constraint or threshold on the time step. 

Finally, the main purpose of this work has been to propose a consistent and specific framework 
for the simulation of polydispersed two-phase flows based on Langevin stochastic equations. To- 
gether with the presentation of the theoretical aspects [5] , it provides a comprehensive description 
of the model and of the numerical ideas. It does not pretend to be the ultimate word in this field 
and much work remains to be done. The Langevin equations still require new developments [8] 
and, on the numerical side, boundary conditions must be properly addressed, see Section [6731 Yet, 
it is hoped that the present framework paves the way for the improvement of current methods as 
well as for the formulation of new ideas. In particular, new hybrid methods may benefit from these 
first steps, by trying to go further into a mixed mean-field/PDF approach within the description 
of the particle phase itself. This will require a good understanding of the consistency between the 
mean field equations satisfied by particle statistical properties and the instantaneous stochastic 
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equations for the trajectories, of their mathematical manipulation and of the issues involved with 
particle/mesh exchange of information. 



Acknowledgements 

Dr. Chibbaro's research was supported by a Marie Curie Transfer of Knowledge fellowship of the 
European Community programme, Contract No. MTKD-CT-2004-509849. The authors would 
like to thank Dr. Mehdi Ouraou for his fruitful collaboration in numerical simulations and Prof. 
Denis Talay for his precious help in developing numerical schemes. 



49 



A Appendix: two-point description 



Here, some additional information is given on the construction of the two-point description, i.e., 
the form of the acceleration term to be added in the Langevin equation describing the velocity 
increments along the trajectory of a fluid particle. Once this is done, it is briefly explained how 
the mean-field (RANS) equations for the fluid can be extracted from this two-point description. 



A.l Model for a two-way coupling term 

In the exact local instantaneous equations for the fluid (the Navier-Stokes equations), a formal 
treatment of the force exerted on the fluid by the discrete particles implies the use of a distribution 
(or density of force) acting on the fluid located in the neighbourhood of the discrete particles in 
order to express the resulting acceleration on nearby fluid particles. This accurate treatment, 
which would result in a multi-point treatment of the discrete phase, is outside the scope of the 
present work. Here, in the frame of the one-point approach, the influence of the discrete particles 
on the fluid is expressed directly in the SDEs, Eqs. (|30|) . with stochastic tools. 

As explained in Section f3. 1.21 for A p ^. s , the underlying force corresponds to the exchange of 
momentum between the fluid and the particles (drag force). The acceleration acting on the fluid 
element surrounding a discrete particle can be obtained as the sum of all elementary accelerations 
(due to the neighbouring particles), i.e., at the discrete particle location x p , the elementary accel- 
eration (Up — U s )/r p is multiplied by \ — (a p p p )/(ctfPf), that is the probable mass of particles 
divided by the probable mass of fluid (since the total force is distributed only on the fluid phase). 

For Ap_>.f , the problem of finding a suitable stochastic model is slightly more difficult since the 
drag force can only be defined in terms of variables attached to the discrete particles (which are 
not defined at the location of a fluid particle). As a consequence, the influence of the neighbouring 
discrete particles on the fluid particle located, at time t, at x = x/(i), is ensured by considering 
that Ap_> / is a random variable given by 

{0 with probability 1 — a p (i,x), 
(124) 
lip with probability a p (t,x), 

where lip is a random variable which plays the role of an ersatz of the Eulerian random variable 
which is formed from the discrete particles at x = x p (t) 

n =-£p Ufl - u p , (125) 

Pi t p 

This random term mimics the reverse force due to the discrete particles and is only non zero when 
the fluid particle is in the close neighbourhood of a discrete particle. In addition, it is required 
that, at x = x/(i), n p and U(t,x) are correlated so that 

r (Hp> = -( Pp / Pf )((u s u p )/Tp), 

\ (HpU) = -(p p /p f )((V s - Up)U s /r p ). 
For further explanations on the modelling of two-way coupling, see Refs. [8) I29j. 



A. 2 Mean-field (RANS) equations for the fluid 

In sample space, Eqs. (f5U| are equivalent (in a weak sense) to a general Fokker-Planck equation for 
the two-point Lagrangian pdf, p r (t; z/, z p ), cf. the correspondence between Eq. (fTt)j) and Eq. (|TT|) . 
It can be shown [8] that the Fokker-Planck equation verified by p r (t;Zf,z p ) is also verified by the 
two-point Eulerian mass density function (mdf ) and therefore by one of its marg 
This mdf is given by F®(t, x; V/) = pj pj (i, x; V /) where pj is the Eulerian distribution function 
of the fluid. The knowledge of the PDE verified by an Eulerian quantity allows us, using classical 
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tools of kinetic theory [SSI US] , to write field equations for the velocity moments of the fluid: the 
PDE verified by Fj is multiplied by a given function of V/, H(Vf). Applying the following 
operator 

a f (t,xi) p f (H(lJ(t,x)) = J n(V f )F k E (t^V f )dV f (127) 

to this PDE gives field equations for any (%(U(i, x)). By replacing H by % = 1, H = Vf,% and 
H = VfjVfj, the continuity equation, the momentum equations and the Reynolds-stress equations 
are obtained, respectively [29]. These equations are given in Table [I] i.e. Eqs. (|128p - (|130|) . 
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m • 


s" 


-3/2 


B* 


approximated value defined by Eq. (|88l) 
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diffusion matrix defined by Eqs. ([79]) and (|12()|) 
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function defined in Table [5] 
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proportionality constant, cf. Eq. (f5T|) 
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proportionality constant, cf. Eq. (|67p 
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drag coefficient 










acceleration defined by Eq. P0|) 
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coefficient for T* L l , cf. Eq. flU) 
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Kolmogorov constant, cf. Eq. (|2"2"|) 
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coefficient defined in Table [4] 


s 






C2 C 


function defined in Tabled] 
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time increment or observation time scale 
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diffusion coefficient or cyclone diameter 
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positive-definite matrix, D = BB T 
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coefficient defined in Table |4] 
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gravitational acceleration 
return-to-equilibrium matrix, cf. Eq. (|23[) 
anisotropy matrix, cf. Eq. (p?5|) 
standard Gaussian random variable, Eq. (|72l) 
standard Gaussian random variable, Eq. (|96|) 
standard Gaussian random variable, Eq. (|73l) 
standard Gaussian random variable, cf. Table |4] 
standard Gaussian random variable, cf. Tabled] 
standard Gaussian random variable, cf. Tabic H] 
stochastic process, cf. Eq. (|G8I) 
function of V/ 

stochastic integral in model problem, cf. Eq 
stochastic integral in model problem, cf. Eq 
turbulent kinetic energy 
modified turbulent kinetic energy 
deposition velocity 

deterministic function, cf. Section [5.2. II 

total number of discrete particles or samples 

number of discrete particles in cell k 

number of discrete particles per cell or per class 

number of discrete particles in cell i 

probability density function (pdf) for Z(t) 

Eulerian fluid distribution function 
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local instantaneous pressure field 
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Reynolds stress tensor 

deterministic function, cf. Section [523] 

source term in Eq. (|129[) 

trace of source term tensor in Eq. (|130[) 

source term in Eq. (|130[) 
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dimensionless time, t + = t (u*) 2 /uf 

characteristic time scale 

fluid Eulerian integral time scale 

fluid seen integral time scale (in Section [SJ 
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fluctuating fluid velocity field 
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local instantaneous fluid velocity field 
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approximated value of U/(i) at t n +i 
velocity of the discrete particles 
approximated value of U p (i) at 
predicted value of U p (f) at t„+i 
approximated value of U p (i) at f n +i 
particle relative velocity 
fluid velocity seen 
approximated value of XJ S (t) at t n 
predicted value of U s (£) at t n+ i 
approximated value of U a (i) at t n +j 
sample space value for U / (t) 
volume of fluid in cell [fc] 
weighting function (continuous form) 
weighting function (discrete form) 
Wiener process 

position in model problem, cf. Eq. 
position of the fluid particles 
position of the discrete particles 
approximated value of x p (i) at t n 
approximated value of x p (i) at i n +i 
stochastic process or deterministic variable, X 
dimensionless distance from wall 
sample space value of Y(t) 
stochastic process or deterministic variable, Y 
set of external variables 
sample space value of Z(t) 
sample space value of Z\(t) 
sample space value of Z*(t) 
state vector 

variable j for fluid particle i 
variable j for discrete particle i 
reduced state vector 
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Greek letters 



a/(t,x) 


volume fraction of fluid 






a p (i,x) 


volume fraction of particles 








constants defined in Eq. (|25|) 
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stochastic process defined by Eq. (|137|) 
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stochastic process defined bv Eq. (|138|) 
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Dirac delta function 








Kronecker's symbol 
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time step 
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characteristic cell size 
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dissipation rate of k 
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energy dissipation for fluid particles 
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Gaussian white noise 
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approximated value of 0£ at i„ 
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dynamic viscosity of fluid 




Pa-s 






kinematic viscosity of fluid 
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standard Gaussian random variable, cf. Eq. (|5ip 








standard Gaussian random variable, cf. Eq. 
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standard Gaussian random variable, cf. Eq. 
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random acceleration defined bv Eq. (|125|l 
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Pf 


density of fluid 
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density of discrete particles 




kg.m" 


-3 


(7 


diffusion matrix in Eq. (|45[) 
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standard deviation of • 
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variance of $ 
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characteristic time scale 
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particle relaxation time 
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dimensionless particle relaxation time 
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approximated value of t p at t n 




s 




~ n+1 
Tp 


predicted value of t p at t n+ i 
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Kolmogorov time scale 
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ratio, x = a s p f /a p p p 










sample space value of $ 










stochastic process defined bv Eq. (|139|) 
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approximated value of Qj(t) at t n 
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Subscripts 

/ continuous phase (fluid) 

p discrete phase (particles) 

s fluid properties sampled along particle trajectories 
Superscripts 



['] 


variable calculated at cell centre i 


[k] 


variable calculated in cell k 


n 


approximated values at t = t n 


n + 1 


approximated values at t = t n + At 


(n), (N) 


variables calculated at particle locations 


r 


reduced information 


T 


transpose of a matrix 


[x] 


variables calculated on the mesh/at cell centres 


+ 


dimensionless quantities 




predicted quantities (numerical schemes) 



Special notation 

{ • } set of variables 

( • ) mathematical expectation 

(•)iv mean value, i.e. (l/JV)S i=1 - 

( • ) a spatial average 

( ■ )n,a approximation of ( • ), spatial average on TV samples 

( • )oo = (-)n with N -> oo, i.e. ( ■ ) 
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( • I ■ ) conditional expectation 

| • | norm of a vector 

d partial derivative 

U bold style for vector notation 

D'/Dt d ■ /dt +(Ui) d ■ /dn 

d ■ (t) time increment, e.g. <fU/(t) = U/(i + dt) — U/(i) 
Abbreviations 



CIC 


Cloud In Cell 


CPU 


Central Processing Unit 


est 


a given constant 


CTE 


Crossing Trajectory Effect 


DNS 


Direct Numerical Simulation 


DSMC 


Direct Simulation Monte Carlo 


LES 


Large Eddy Simulation 


NGP 


Nearest Grid Point 


ODE 


Ordinary Differential Equation 


PDE 


Partial Differential Equation 


pdf/PDF 


Probability Density Function 


PIC 


Particle In Cell 


RANS 


Reynolds-Averaged Navier-Stokes 


RHS 


Right-Hand Side 


r.m.s. 


root-mean square 


RSM 


Reynolds Stress Models 


SDE 


Stochastic Differential Equation 


SPH 


Smoothed Particle Hydrodynamics 


VRT 


Variance Reduction Technique 
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Figure 1: Weighting functions of different orders used for averaging: (a) top hat (constant) or 
Nearest-Grid-Point (NGP); (b) linear or Cloud-In-Cell (CIC); (c) piecewise quadratic. 
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Figure 2: CIC scheme used to compute cell averages for (a) particle density and (b) a linear function 
attached to particles (the linear function is defined by its slope A$ and <&o which is the value of $ 
at x = ajt 8+1 l). The cells are characterised by the cell boundary coordinates Xq 1 , Xq , Xq 1 , . . . 
and by the coordinates of the cell centres x^ _1 \x^ ,x^ l+1 \ .... Axt and Axi+i represent the cell 
sizes. 
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Figure 3: Computations of the r.m.s. of particle density on a mesh using NGP and CIC averaging 
: (a) uniform mesh, (b) non-uniform mesh. NGP method (□), CIC method (■). 
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Figure 4: The r.m.s. value of the mean of a linear deterministic function computed on a mesh 
using NGP and CIC averaging : (a) uniform mesh, (b) non- uniform mesh. NGP method (□), CIC 
method (■). 
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Figure 5: Normalised r.m.s. of a linear random function computed on a mesh using NGP and CIC 
averaging: (a) uniform mesh, (b) nonuniform mesh. NGP method (□), CIC method (■). 
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Figure 6: Cells and particles - schematic plot. Mirror particles (□), corresponding to (■), are 
added outside of the computational domain. The values of variables attached to them (like veloc- 
ity) correspond to those of their "host" particles in border cells. Dashed lines delimit cells and 
dotted lines indicate where mirror particles are needed. 
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Figure 7: Temporal evolution of the mean (lower plots), (<£>), and variance (upper plots), (o$), 
profiles for quadratic initial mean profile and constant initial non-zero variance. Left plots: CIC; 
right plots: NGP. Three successive time instants (solid, dotted and dashed lines). 
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Figure 8: Temporal evolution of the r.m.s. of $ in scaled time; NGP computation over 1000 
cells, (a) Results for N pc = 20, 25, 40, 80, 160, and 320, respectively (from the lowest to the highest 
curve); (b) the r.m.s. of $ with its standard deviation for N pc = 25 (dotted line with o) and 
N pc = 320 (solid line with •). 
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Figure 9: (a) Analysis of the bias of the r.m.s. of <J> as a function of the number of particles per 
cell N pc for three different NGP computations; dotted line: -1 slope, (b) PDF of the r.m.s. of $ 
at t/T = 3 in the NGP computation over 5000 cells. Solid line: iVp G = 20, dashed line: N pc = 70. 
Results are smoothed histograms (dashed lines). 
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Figure 10: Definition of the geometry of a cyclone of the Stairmand design. All dimensions are 
given as a function of the diameter D of the cyclone. 
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Figure 11: Gas mean axial velocity profiles at two different heights: z = 0.36 m (left) and 
z = 0.57m (right). The velocities are plotted as functions of the radius R. Continuous lines: 
computations. Experimental data: A. 
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Z = 0.36 m Z = 0.57 m 




Figure 12: Gas mean tangential velocity profiles at two different heights: z = 0.36 m (left) and 
z = 0.57m (right). The velocities are plotted as functions of the radius R. Continuous lines: 
computations. Experimental data: A. 



72 



order 2, At = 10" 4 order 2, At = 1 




/ / 
2 0.2 - / 




d p (jim) 



order 1 , (At x 1 







P*"' At=1 -- 




// At 0.5 




1 At 0.1 







dp (urn) 




dp (jim) 



comparison order 1 /order 2 



order 2, At=10" 
order 1, At=5.10"' 




Figure 13: Sensitivity analysis for the second order and first order schemes. The sensitivity analysis 
is carried out for the number of particles per class, N pc , and the time step, At (with N pc — 400). 
At last, comparison of the time steps for the first and second order schemes for identical numerical 
results (bottom right corner). All numerical results represent the selectivity curve which gives the 



efficiency of the cyclone, e as a function of the particle diameter, d 
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x (gravity) 

Figure 14: The 'Hercule' experimental setup. The mean streamlines are shown for the fluid (solid 
lines) and the particles (dashed lines). Two stagnation points in the fluid flow can be observed 
(Si and S2). Experimental data is available for radial profiles of different statistical quantities at 
five axial distances downstream of the injection (x = 0.08,0.16,0.24,0.32,0.40m) (experimental 
data is also available on the symmetry axis). 



74 



x = 0.08 m x = 0.16 m 




0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.02 0.04 0.06 0.08 0.1 0.12 0.14 

r (m) r (m) 



Figure 15: Radial profiles of the particle axial velocity at x = 0.08 m and x = 0.16 m. Mean 
velocities (top) and fluctuating velocities (bottom). Continuous lines: computations (first order 
scheme). Experimental data: ■. 
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Figure 16: Radial profiles of the particle radial velocity at x = 0.08 m and x = 0.16 m. Mean 
velocities (top) and fluctuating velocities (bottom). Continuous lines: computations (first order 
scheme). Experimental data: ■. 
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Figure 17: Radial profiles of the fluid axial velocity at x 




0.02 0.04 0.06 0.08 0.1 0.12 0.14 
z (m) 

0.08 m and x = 0.16 m (mean velocities 



only). Continuous lines: computations (first order scheme). Experimental data: 
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Figure 18: Left-hand side: snapshot of the flow, at a given time, where the particles are coloured 
by their residence time in the computational domain. Right-hand side (upper corner): probability 
density functions of the particle residence time inside the domain (close to the injection and close 
to the outlet of the domain) . Right-hand side (bottom corner) : axial velocity (horizontal axis) of 
the particles (coloured by their residence time in the computational domain) close to the inlet at 
different time steps (vertical axis). 78 
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Figure 19: Sensitivity analysis : deposition velocity computed with different time steps At and 
with the first and second order schemes. At = 10~ 4 s (o) and At = 10~ 5 s (o) with the first order 
scheme. At = 10 _4 s (□) and At = 10 -5 s (A) with the second order scheme. 
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Figure 20: Sensitivity analysis : deposition velocity computed with different number of particles 
per class N pc for the same fluid mean-fields (N pc = 500 (o), N pc = 1000 (□) and N pc = 5000 (o)). 
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Figure 21: Deposition velocity with different mean fields for the fluid. Fluid mean-fields calculated 
with a k — e model (o). Fluid mean-fields computed with a Rtj — e model (□). Experimental data 

(v). 
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Figure 22: Deposition velocity with different mean fields for the fluid. First, mean field obtained 
by computation with a standard k — e model (o). Second, mean field where (Ui) is computed 
with the law-of-the-wall equations and where the values k and e are that of the computations (o). 
Third, (Ui) is still given by the law-of-the-wall, and k and e are curve-fitted to DNS data [65] (A). 
Experimental data (V). 
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Table 1: Complete mean-field (RANS)/PDF model. 



Mean-field (RANS) equations for the fluid. 
Continuity equation: 

_(a /P/ ) = -a /P/ -^ (128) 
Momentum equation: 

?rM) = --^P- —^(afPfiuiU^+xW^i - U s ,, t )/r p ) (129) 
Dt p f oxi a.f pf oxj 

Reynolds stress equation: 

D Id, , ,. . .d(Uj) , ,d(Ui) 

Di {UiUj) = -^Tfdx-^ Pf{UlU > Uk)) ~ M ^xV ^ Uk) ^xV 

+ G jk (uiU k ) + G. lk (ujU k ) + C (e)6ij 

+ X(— [(U p ,i ~ U s>i )U S!j + {U p ,j - U s ,j)U s ,i\) (130) 

Tp 

D d d a„pp 

with —— = — + \Uk) — and x = - J ^ L 

Dt at oxk a f Pf 

SDEs for the discrete particles. 

dx p j(t) = U Pt i dt, 
dUp At) = — ^ — — dt + gi dt, 

Tp 

dU s>i (t) = A Sil dt + Ap^ s>i dt + B Stij dWj(t), 

Pf uxi uxj ± Li 

B 2 s>i = (e) (Cohk/k+^ibik/k-l) 

-Ap— > s ,i = \{Us,i Up,i)/Tp 
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Table 2: Analytical solutions to system (QUI) for time-independent coefficients. 



x P ,i(t) = x P!i (to) + U p ,i(ta)T p [l - exp(-At/r p )] + U s ,i(t ) 9 l {T l [l - exp(-At/T,)] 
+ T p [exp(-Ai/r p ) - 1]} + [Q T t ]{At - r p [l - exp(-At/r p )] 

- Oi(Ti[l - exp(-At/T,)] + T p [exp(-At/r p ) - 1])} + O^t) (131) 
with 6i =Ti/ (Ti - t p ) 

U Pii (t) = U Pii (to) exp(-At/r p ) + U %l (t ) 9 t [exp(— At/Tj) - cxp(-At/r p )] 
+ [C t Ti\{[l ~ exp(-At/T p )) - 9i[exp(-At/Ti) - cxp(-At/r p )]} 
+ T l (t) (132) 

U s ,i(t) = C/ S)i (i ) exp(-At/T J ) + Q T % [\ - exp(-Ai/T i )] + 7i (t) (133) 

The stochastic integrals 7i(t), r,(i), are given by: 



7i (t) = B< exp(-t/T<) / exp(s/Ti)dWi(s), (134) 

J to 

1 /■< 

Ti(t) = —exp(-t/T p ) exp(s/T p )%(s)ds, (135) 

T P -'to 

fii(t) = / ri(s)ds. (136) 
•/to 

By resorting to stochastic integration by parts, 7i(t), £li(t) can be written: 

7i(t) = Bi exp(—t/Ti)Ii t i, (137) 
(t) = 0, [exp(-t/Ti) Ji,< - exp(-t/r p ) J a ,i] , (138) 

n i (t) = fl < 5 1 -{(ii-T J ,)j 3 ,< 

- [r i exp(-VTi)7i,i-r p exp(-t/r p )/2,i]}, (139) 
with J M = / cxp(s/Ti)dWi(s), J 2)i = / exp(s/T p )dWi(s) 

J t J t 



and I 3<i = / dW*(s). 
Jt 
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Tabic 3: Derivation of the covariance matrix for constant coefficients. 



(jf(t)) = Bf ^ [1 - exp(-2At/T 2 )] where &* = B\ (140) 

(T?(t)} = 0? - exp(-2AVT i )] - ^-[1 - exp(-Ai/T,) exp(-At/r p )] 

+ ^[l-exp(-2AVr p )]} (141) 

(0 2 (t)) - (Tj - r p ) 2 At + ?[1 - exp(-2At/T i )] + ^[1 - exp(-2At/r J) )] 



S| (9? x * w/ v P ' 2 ^ v ' i/J 2 

- 2I?(T< - Tp )[l - exp(-AVIi)] + 2r 2 (T l - r p )[l - cxp(-At/r p )] 

T 2 r 2 

- [1 - cx P (- At/Ti) exp(- At/rp)] (142) 



T, 



( 7< (t) r<(t)) - 5? 0,2} f i[l - exp(-2Ai/2})] - jjrj— [1 - exp(-At/2}) exp(-At/T p )]\ (143) 

12 J-i+T p j 

( 7i (t) Jli(t)) = B 2 4 2} {(2} - r p )[l - exp(-Ai/2})] - y [1 - exp(-2At/2})] 
+ TfTT — P- - cxp{- At/Ti) cxp(-At/T p )} \ (144) 



2} + r p 



(^(i) fii(t)) = (2} - r p ){2}[l - exp(-Af/2})] - Tp [l ~ exp(-At/Tp)]} 



2 fl 2 



2 [1 - exp(-2Ai/2})] - -j- [1 - cxp(-2At/rp)] 
+ 2}r p [1 - exp(-At/2}) exp(-At/r p )] (145) 
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Tabic 4: Weak first-order scheme (Eulcr scheme) 



Numerical integration of the system: 



= 17£ e^-Af/I?) + K"Cf][l - exp(-Ai/Tf )] + 7 f , 



The coefficients A x , i?i , Ci , D\ and £1 are given by: 



A 1= r;[l-exp(-AVr;)], 

B X = 9? [T? (1 - expC-At/TI 1 ) - Ai] with 0? = Tf/(Tf - r p "), 
Ci = Ai-Ai-Bi, 

D x = 0?[exp(-Ai/Tf ) - exp(-At/r;)], 
E x = l-cxp(-At/r; 1 ). 



The stochastic integrals 7™, f2", T" are simulated by: 



where Gi^, G2,i, Gz,i are independent 7V(0, 1) random variables. 
The coefficients Pn, P 2 i, -P22, P31, P32, P33 are defined as: 



P3i= 7==? > P*2 = 1 ^mr?)-p 21 p 31 ), P33 = v / ((rr) 2 )-^ 2 i-^3 2 2 )- 



7? = ^!^, 

i? - p 3 i ft,* + -P32 £ 2 . 4 + P33 G 3 , 



Pa 
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Tabic 5: Weak second-order scheme 



Prediction step: Euler scheme, see Table |U 



Correction step: 



U, 



n+l 



1 



XJli exp(-At/7£) + - V%. exp(-AVf; +i ) 



1 



1 



~n+l 

p 



Ti ) 



+ A 2c (t£, Tf) [T?C?] + B 2c (t£ +1 , f t n+1 ) [f« +1 C? +1 ] 

~ n+l\r ~ n+l An+1] , fm+1 



If. 



n+l 

S,2 



+ ^ 2 (At,r;)[r;^ 1 ] +S 2 (Ai,f p " +1 )[T p " + M™ +1 ] +f- 
+ B 2 (At, f t n+1 ) [f/ l+1 C™+ 1 ] + 7f +1 . 



i 0£ exp(-At/lT) + \ Ki cx P (-Ai/f; n+1 ) + A 2 (Ai, 1?) pTC?] 



The coefficients A 2 , B 2 , A 2c , i? 2c et C 2c are defined as: 
A 2 (At,x) = - exp(-Ai/x) + [1 - exp(-Af/x)][Ai/a:], 
B 2 (At,x) = 1 - [1 -exp(-Ai/x)][Ai/a;], 

A 2c (x, y) = - eoqp(-At/s) + [(a; + y)/At] [1 - exp(-Ai/x)] - (1 + y/M) C 2c (x, y), 
B 2c (x, y) = l — [(x + y)/At][l - exp(-At/x)] + (y/At) C 2c (x, y), 
C 2c {x,y) = \y/(y-x)][exp(-£t/y)-exp(-At/x)]. 

The stochastic integrals 7™ +1 and arc simulated as follows: 



with 



n+l 



-[l-cx P (-2AV^ +i )] g lti , 
V z 

1 - exp(-2 At/f? +1 )\ B* = A 2 (2 At, T/ l+1 ) yJ(B?) 2 + 

s 2 (2Af,r» +1 ) v /(^r +1 ) 2 - 



(T™ +1; v ,l+1 \ / /p»+Un+l\2 



((7r +1 ) 2 > 



((7r +1 ) 2 > 



with (rr i 7r i > = <r^X^,Tr +1 , J B;) and ((r™ ) ) = (r^)(r™, T" , B*) 
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Table 6: Mean near- wall (y + < 30) residence time for different diameters (the simulation is 
carried out with the exact frozen field). The residence time is given in non-dimensional form (it 
is normalised with the viscous time scale, vf /(u*) 2 , that is t + = t (u*) 2 jvf). 



T+ 


dp (fim) 


t + (wall units) 


0.2 


1.4 


29.5 


0.4 


2.0 


29.9 


0.9 


2.9 


28.7 


1.9 


4.3 


30.9 


3.5 


5.8 


31.5 


6.4 


7.8 


31.6 


13.2 


11.2 


35.4 


29.6 


16.8 


40.6 


122.7 


34.2 


55.3 


492.2 


68.5 


96.8 
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